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BACKGROUND OF THE INVENTION 

a . Field of Invention 

The invention relates generally to computerized 
15 techniques for processing data obtained from radar to track 
multiple discrete objects. 

b . Description of the Background 

There are many situations where the courses of multiple 
objects in a region must be tracked. Typically, radar is used 

2 0 to scan the region and generate discrete images or "snapshots" 
based on sets of returns or observations. In some types of 
tracking systems, all the returns from any one object are 
represented in an image as a single point unrelated to the 
shape or size of the objects. "Tracking" is the process of 

25 identifying a sequence of points from a respective sequence of 
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the images that represents the motion of an object. The 
tracking problem is difficult when there are multiple closely 
spaced objects because the objects can change speed and 
direction rapidly and move into and out of the line of sight 
5 for other objects. The problem is exacerbated because each 
set of returns may result from noise as well as echoes from 
the actual objects. The returns resulting from the noise are 
also called false positives. Likewise, the radar will not 
detect all echoes from the actual objects and this phenomena 

10 is called a false negative or "missed detect" error. For 
tracking airborne objects, a large distance between the radar 
and the objects diminishes the signal to noise ratio so the 
number of false positives and false negatives can be high. 
For robotic applications, the power of the radar is low and as 

15 a result, the signal to noise ratio can also be low and the 
number of false positives and false negatives high. 

In view of the proximity of the objects to one another, 
varied motion of the objects and false positives and false 
negatives, multiple sequential images should be analyzed 

20 collectively to obtain enough information to properly assign 
the points to the proper tracks. Naturally, the larger the 
number of images that are analyzed, the greater the amount of 
information that must be processed. 
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While identifying the track of an object, a kinematic 
model may be generated describing the location, velocity and 
acceleration of the object. Such a model provides the means 
by which the future motion of the object can be predicted. 
5 Based upon such a prediction, appropriate action may be 
initiated. For example, in a military application there is a 
need to track multiple enemy aircraft or missiles in a region 
to predict their objective, plan responses and intercept them. 
Alternatively, in a commercial air traffic control application 

10 there is a need to track multiple commercial aircraft around 
an airport to predict their future courses and avoid 
collision. Further, these and other applications, such as 
robotic applications, may use radar, sonar, infrared or other 
object detecting radiation bandwidths for tracking objects. 

15 In particular, in robotic applications reflected radiation can 
be used to track a single object which moves relative to the 
robot (or vice versa) so the robot can work on the object. 

Consider the very simple example of two objects being 
tracked and no false positives or false negatives. The radar, 

20 after scanning at time t lf reports objects at two locations in 
a first observation set. That is, the radar returns a set of 
two observations {o llf o 12 ] . At time t 2 the radar returns a 
similar set of two observations {o 21/ o 22 ) from a second 
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observation set. Suppose from prior processing that track 
data for two tracks T 2 and T 2 includes the locations at t 0 of 
two objects. Track T x may be extended through the points in 
the two sets of observations in any of four ways, as may track 
5 T 2 . The possible extensions of T 2 can be described as: {T x , 
o llf o 21 } , {T lf o llf o 22 ) , {T lf o 12f o 21 ) and {T lf o 12 , o 22 ) . Tracks 
can likewise be extended from T 2 in four possible ways, 
including {T 2 , o 12/ o 21 } . Figure 1 illustrates these five (out 
of eight) possible tracks (to simplify the problem for 

10 purposes of explanation) . The five track extensions are 
labeled h llf h 12 , h 13f h 14/ and h 21 wherein h ±1 is derived from [T lf 
°iif °2i] / &12 i s derived from {T 2/ o llf o 22 ] , h 13 is derived from 
{T 2/ o 12 , o 21 } , h 14 is derived from [T lf o 12/ o 22 ] , and h 21 is 
derived from {T 2r o llf o 21 }. The problem of determining which 

15 such track extensions are the most likely or optimal is 
hereinafter known as the assignment problem. 

It is known from prior art to determine a figure of merit 
or cost for assigning each of the points in the images to a 
track. The figure of merit or cost is based on the likelihood 

2 0 that the point is actually part of the track. For example, 
the figure of merit or cost may be based on the distance from 
the point to an extrapolation of the track. Figure 1 
illustrates costs 5 21 and 5 22 for hypothetical extension h 21 and 
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modeled target characteristics. The function to calculate the 
cost will normally incorporate detailed characteristics of the 
sensor, such as probability of measurement error, and track 
characteristics, such as likelihood of track maneuver. 
5 Figure 2 illustrates a two by two by two matrix, c, that 

contains the costs for each point in relation to each possible 
track. The cost matrix is indexed along one axis by the track 
number, along another axis by the image number and along the 
third axis by a point number. Thus, each position in the cost 

10 matrix lists the cost for a unique combination of points and 
a track, one point from each image. Figure 2 also illustrates 
a {0, 1} assignment matrix, z, which is defined with the same 
dimensions as the cost matrix. Setting a position in the 
assignment matrix to "one" means that the equivalent position 

15 in the cost matrix is selected into the solution. The 
illustrated solution matrix selects the {h 14 , h 21 ) solution 
previously described. Note that for the above example of two 
tracks and two snapshots, the resulting cost and assignment 
matrices are three-dimensional. As used in this patent 

20 application, the term "dimension" means the number of axes in 
the cost or assignment matrix while size refers to the number 
of elements along a typical axis. The costs and assignments 
have been grouped in matrices to facilitate computation. 
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A solution to the assignment problem satisfies two 
constraints - first, the sum of the associated costs for 
assigning points to a track extension is minimized and, 
second, if no false positives or false negatives exist, then 
5 each point is assigned to one and only one track. 

When false positives exist, however, additional 
hypothetical track extensions incorporating the false 
positives will be generated. Further note that the random 
locations of false positives will, in general, not fit well 

10 with true data and such additional hypothetical track 
extensions will result in higher costs. Also note that when 
false negative errors exist, then the size of the cost matrix 
must grow to include hypothetical track extensions formulated 
with "gaps" (i.e., data omissions where there should be 

15 legitimate observation data) for the false negatives. Thus, 
the second criteria must be weakened to reflect false 
positives not being as signed and also to permit the gap 
filler to be multiply assigned. With hypothetical cost 
calculated in this manner then the foregoing criteria for 

20 minimization will tend to materialize the false negatives and 
avoid the false positives. 

For a three-dimensional problem, as is illustrated in 
Figure 1, but with N x (initial) tracks, N 2 observations in scan 
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1, N 3 observations in scan 2, false positives and negatives 
assumed, the assignment problem can be formulated as: 

W- n 0 w, 



(a) Minimize: 2J zJ X/ c i i i z ± ± ± 

i^o 1^0 i 3 =o ll12 ^ 11X2 3 



(b) Subject to: 



i 2 =i i 3 -i ^3 



5 (c) }2 22 i 2 = 1,...N 2 , 

i 3 =i 1 2 3 

VI i 2 =l 

(e) z ii± € {0,1} V z nj , 



10 where "c" is the cost and "z" is a point or observation 
assignment, as in Figure 2. 

The minimization equation or equivalently objective 
function (0.1) (a) specifies the sum of the element by element 
product of the c and z matrices. The summation includes 

15 hypothesis representations z ili2i3 with observation number zero 
being the gap filler observation. Equation (0.1) (b) requires 
that each of the tracks T 2 , . . . , T N be extended by one and only 
one hypothesis. Equation (0.1) (c) relates to each point or 
observation in the first observation set and requires that 

2 0 each such observation, except the gap filler, can only 
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associate with one track but, because of the "less than" 
condition, it might not associate with any track. Equation 
(0.1) (d) is like (0.1) (c) except that it is applicable to the 
second observation set. Equation (0.1) (e) requires that 
5 elements of the solution matrix z be limited to the zero and 
one values. 

The only known method to solve Problem Formulation (0.1) 
exactly is a method called "Branch and Bound." This method 
provides a systematic ordering of the potential solutions so 

10 that solutions with a same partial solution are accessible via 
a branch of a tree describing all possible solutions whereby 
the cost of unexamined solutions on a branch are incrementally 
developed as the cost for other solutions on the branch are 
determined. When the developing cost grows to exceed the 

15 previously known minimal cost (i.e., the bound) then 
enumeration of the tree branch terminates. Evaluation 
continues with a new branch. If evaluation of the cost of a 
particular branch completes, then that branch has lower cost 
than the previous bound so the new cost replaces the old 

2 0 bound. When all possible branches are evaluated or eliminated 
then the branch that had resulted in the last used bound is 
the solution. If we assume that N ± = N 2 = N 3 = n and that 
branches typically evaluate to half there full length, then 
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workload associated with "branch and bound" is proportional to 
(n!|-^!) 2 . This workload is unsuited to real time evaluation. 

The Branch and Bound algorithm has been used in past 
research on the Traveling Salesman Problem. Messrs. Held and 
5 Karp showed that if the set of constraints was relaxed by a 
method of Lagrangian Multipliers (described in more detail 
below) then tight lower bounds could be developed in advance 
of enumerating any branch of the potential solution. By 
initiating the branch and bound algorithm with such a tight 
10 lower bound, significant performance improvements result in 
that branches will typically evaluate to less than half their 
full length. 

Messrs. Frieze and Yadagar in dealing with a problem 
related to scheduling combinations of resources, as in job, 

15 worker and work site, showed that Problem Formulation (0.1) 
applied. They further described a solution method based upon 
an extension of the Lagrangian Relaxation previously 
mentioned. The two critical extensions provided by Messrs. 
Frieze and Yadagar were: (1) an iterative procedure that 

2 0 permitted the lower bound on the solution to be improved (by 
"hill climbing" described below) and (2) the recognition that 
when the lower bound of the relaxed problem was maximized, 
then there existed a method to recover the solution of the 
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non-relaxed problem in most cases using parameters determined 
at the maximum. The procedures attributed to Messrs. Frieze 
and Yadagar are only applicable to the three-dimensional 
problem posed by Problem Formulation (0.1) and where the cost 
matrix is fully populated. However, tracking multiple 
airborne objects usually requires solution of a much higher 
dimensional problem. 

Figures 1 and 2 illustrate an example where "look ahead" 
data from the second image improved the assignment accuracy 
for the first image. Without the look ahead, and based only 
upon a simple nearest neighbor approach, the assignments in 
the first set would have been reversed. Problem Formulation 
(0.1) and the prior art only permit looking ahead one image. 
In the prior art it was known that the accuracy of assignments 
will improve if the process looks further ahead, however no 
practical method to optimally incorporate look ahead data 
existed. Many real radar tracking problems involve hundreds 
of tracks, thousands of observations per observation set and 
matrices with dimensions in the range of 3 to 25 including 
many images of look ahead. 

It was also known that the data assignment problem may be 
simplified (without reducing the dimension of the assignment 
problem) by eliminating from consideration for each track 
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those points which, after considering estimated limits of 
speed and turning ability of the objects, could not physically 
be part of the track. One such technique, denoted hereinafter 
the "cone method, 11 defines a cone as a continuation of each 
5 previously determined track with the apex of the cone at the 
end of the previously defined track. The length of the cone 
is based on the estimated maximum speed of the object and the 
size of the arc of the cone is based on the estimated maximum 
turning ability of the object. Thus, the cone defines a 

10 region outside of which no point could physically be part of 
the respective track. For any such points outside of the 
cones, an infinite number could be put in the cost matrix and 
a zero could be preassigned in the assignment matrix. it was 
known for the tracking problem that these elements will be 

15 very common in the cost and selection matrices (so these 
matrices are "sparse"). 

It was also known in the prior art that one or more 
tracks which are substantially separated geographically from 
the other tracks can be separated also in the assignment 

2 0 problem. This is done by examining the distances from each 
point to the various possible tracks. If the distances from 
one set of points are reasonably short only in relation to one 
track, then they are assigned to that track and not further 
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considered with the remainder of the points. Similarly, if a 
larger group of points can only be assigned to a few tracks, 
then the group is considered a different assignment problem. 
Because the complexity of assignment problems increases 
5 dramatically with the number of possible tracks and the total 
number of points in each matrix, this partitioning of the 
group of points into a separate assignment problem and removal 
of these points from the matrices for the remaining points, 
substantially reduces the complexity of the overall assignment 
10 problem. 

A previously known Multiple Hypothesis Testing (MHT) 
algorithm (see Chapter 10 of S.S. Blackman. Multiple-Target 
Tracking with Radar Applications. Artech House, Norwood, MA, 
1986.) related to formulation and scoring. The MHT procedure 

15 describes how to formulate the sparse set of all reasonable 
extension hypothesis (for Figure 1 the set {h 1± . . .h 24 } ) and how 
to calculate a cost of the hypothesis {r l7 o ljf o 2k } based upon 
the previously calculated cost for hypothesis {T u o 2j } . The 
experience with the MHT algorithm, known in the prior art, is 

20 the basis for the assertion that look ahead through k sets of 
observations results in improved assignment of observations 
from the first set to the track. 
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In theory, the MHT procedure uses the extendable costing 
procedure to defer assignment decision until the accumulated 
evidence supporting the assignment becomes overwhelming. When 
it makes the assignment decision, it then eliminates all 
5 potential assignments invalidated by the decision in a process 
called "pruning the tree." In practice, the MHT hypothesis 
tree is limited to a fixed number of generations and the 
overwhelming evidence rule is replaced by a most likely and 
feasible rule. This rule considers each track independently 

10 of others and is therefore a local decision rule. 

A general object of the present invention is to provide 
an efficient and accurate process for assigning each point 
object in a region from multiple images to a proper track and 
then taking an action based upon the assignments. 

15 A more specific object of the present invention is to 

provide a technique of the foregoing type which determines the 
solution of a k~ dimensional assignment problem where "Jc" is 
greater than or equal to three. 

2 0 SUMMARY OF THE INVENTION 

The present invention relates to a method and apparatus 

for tracking objects. In particular, the present invention 

tracks movement or traj ectories of obj ects by analyzing 
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radiation reflected from the objects, the invention being 
especially useful for real-time tracking in noisy 
environments . 

In providing such a tracking capability, a region 
5 containing the objects is repeatedly scanned to generate a 
multiplicity of sequential images or data observation sets of 
the region. One or more points (or equivalently 

observations) , in each of the images or observation sets are 
detected wherein each such observation either corresponds to 

10 an actual location of an object or is an erroneous data point 
due to noise. Subsequently, for each observation detected, 
figures of merit or costs are determined for assigning the 
observation to each of a plurality of previously determined 
tracks. * Afterwards, a first optimization problem is 

15 specified which includes: 

(a) a first objective function for relating the above 
mentioned costs to potential track extensions through the 
detected observations (or simply observations) ; and 

(b) a first collection of constraint sets wherein each 
2 0 constraint set includes constraints to be satisfied by 

the observations in a particular scan to which the 
constraint set is related. In general, there is a 
constraint for each observation of the scan, wherein the 
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constraint indicates the number of track extensions to 
which the observation may belong. 

In particular, the first optimization problem is 
formulated, generated or defined as an M-dimensional 
assignment problem wherein there are M constraint sets in the 
first collection of constraint sets (i.e., there are M scans 
being examined) and the first objective function minimizes a 
total cost for assigning observations to various track 
extensions wherein terms are included in the cost, such that 
the terms have the figures of merit or costs for hypothesized 
combinations of assignments of the observations to the 
tracks. Subsequently, the formulated M-dimensional assignment 
problem is solved by reducing the complexity of the problem by 
generating one or more optimization problems each having a 
lower dimension and then solving each lower dimension 
optimization problem. That is, the M-dimensional assignment 
problem is solved by solving a plurality of optimization 
problems each having a lower number of constraint sets. 

The reduction of the M-dimensional assignment problem to 
a lower dimensioned problem is accomplished by relaxing the 
constraints on the points of one or more scans thereby 
permitting these points to be assigned to more than one track 
extension. In relaxing the constraints, terms having penalty 
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factors are added into the objective function thereby 
increasing the total cost of an assignment when one or more 
points are assigned to more than one track. Thus, the 
reduction in complexity by this relaxation process is 
5 iteratively repeated until a sufficiently low dimension is 
attained such that the lower dimensional problem may be solved 
directly by known techniques. 

In one embodiment of the invention, each ^-dimensional 
assignment problem (2 < k < M) is iteratively reduced to a (k- 

10 1)- dimensional problem until a two-dimensional problem is 
specified or formulated. Subsequently, the two-dimensional 
problem formulated is solved directly and a "recovery" 
technique is used to iteratively recover an optimal or near- 
optimal solution to each Jc-dimensional problem from a derived 

15 (A:- 1) -dimensional problem k= 2,3,4,...M. 

In performing each recovery step (of obtaining a solution 
to a Jc-dimensional problem using a solution to a (Jc-1) - 
dimensional problem) , an auxiliary function is utilized. In 
particular, to recover an optimal or near-optimal solution to 

20 a ic-dimensional problem, an auxiliary function, W k . x 
(7c=4,5, . . . ,M) , is specified and a region or domain is 
determined wherein this function is maximized, whereby values 
of the region determine the penalty factors of the (Jc-1) - 
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dimensional problem such that another two-dimensional problem 
can be formulated which determines a solution to the k- 
dimensional problem using the penalty factors of the {k-1) - 
dimensional problem. 
5 Each Auxiliary function ¥ kf Jc=3 , . . . ,AT-1, is a function 

of both lower dimensional problems, penalty factors, and a 
solution at the dimension k at which the penalized cost 
function is solved directly (typically a two-dimensional 
problem) . Further, in determining, for auxiliary function W kf 

10 a peak region, a gradient of the auxiliary function is 
determined, and utilized to identify the peak region. Thus, 
gradients are used for each of the approximation functions ¥ 3 , 
¥ 4f . . . , ¥ M _ X which are sequentially formulated and used in 
determining penalty factors until M-l is used in determining 

15 the penalty factors for the (M-l) -dimensional problem. 
Subsequently, once the M-dimensional problem is solved (using 
a two-dimensional problem to go from an (M-l) -dimensional 
solution to an M-dimensional solution) , one or more of the 
following actions are taken based on the track assignments: 

2 0 sending a warning to aircraft or a ground or sea facility, 
controlling air traffic, controlling anti-aircraft or anti- 
missile equipment, taking evasive action, working on one of 
the objects. 
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According to one feature of this first embodiment of the 
present invention, the following steps are also performed 
before the step of defining the auxiliary function. A 
preliminary auxiliary function is defined for each of the 
5 lower dimensional problems having a dimension equal or one 
greater than the dimension at which the penalized cost 
function is solved directly. The preliminary auxiliary 
function is a function of lower order penalty values and a 
solution at the dimension at which the penalized cost function 

10 was solved directly. In determining a gradient of the 
preliminary auxiliary function, step in the direction of the 
gradient to identify a peak region of the preliminary 
auxiliary function and determine penalty factors at the peak 
region. Iteratively repeat the defining, gradient 

15 determining, stepping and peak determining steps to define 
auxiliary functions at successively higher dimensions until 
the auxiliary function dimension (Jc-1) is determined. In an 
alternative second embodiment of the present invention, 
instead of reducing the dimentionality of the M-dimensional 

2 0 assignment problem by a single dimension at a time, a 
plurality of dimensions are relaxed simultaneously. This new 
strategy has the advantage that when the M-dimensional problem 
is relaxed directly to a two-dimensional assignment problem, 
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then all computations may be performed precisely without 
utilizing an auxiliary function such as W k as in the first 
embodiment. More particularly, the second embodiment solves 
the first optimization problem (i.e., the M-dimensional 
5 assignment problem) by specifying (i.e., creating, generating, 
formulating and/or defining) a second optimization problem. 
The second optimization problem includes a second objective 
function and a second collection of constraint sets wherein: 

a) the second objective function is a combination of the 
10 first objective function and penalty factors or terms 

determined for incorporating (M-m) -constraint sets of the 
first optimization problem into the second objective 
function; 

b) the constraint sets of the second collection include only 
15 m of the constraint sets of the first collection of 

constraints, 

wherein 2 < m < M -1. Note that, once the second optimization 
problem has been specified or formulated, an optimal or near- 
optimal solution is determined and that solution is used in 
20 specifying (i.e, creating, generating, formulating and/or 
defining) a third optimization problem of {M-m) dimensions (or 
equivalently constraint sets) . The third optimization problem 
is subsequently solved by decomposing it using the same 
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procedure of this second embodiment as was used to decompose 
the first optimization problem above. Thus, a plurality of 
instantiations of the third optimization problem are 
specified, each successive instantiation having a lower number 
5 of dimensions, until an instance of the third optimization 
problem is a two-dimensional assignment problem which can be 
solved directly. Subsequently, whenever an instance of the 
third optimization problem is solved, the solution is used to 
recover a solution to the instance of the first optimization 

10 problem from which this instance of the third optimization was 
derived. Thus , an optimal or near-optimal solution to the 
original first optimization problem may be recovered through 
iteration of the above steps. 

As mentioned above, the second embodiment of the present 

15 invention is especially advantageous when m =2, since in this 
case all computations may be performed precisely and without 
utilizing auxiliary functions. 

Other features and benefits of the present invention will 
become apparent from the detailed description with the 

2 0 accompanying drawings contained hereinafter. 
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BRIEF DESCRIPTION OF THE DRAWINGS 
Figure 1 is a graph of images or data sets generated by 
a scan of a region and possible tracks within the images or 
data sets according to the prior art. 
5 Figure 2 illustrates cost and assignment matrices for the 

data sets of Figure 1 according to the prior art. 

Figure 3 is a block diagram of the present invention. 
Figure 4 is a flow chart of a process according to the 
prior art for solving a three-dimensional assignment problem. 
10 Figure 5 is a flow chart of a process according to the 

present invention for solving a Jc-dimensional assignment 
problem where n k n is greater than or equal to 3. 

Figure 6 is a graph of various functions used to explain 
the present invention. 
15 Figure 7 is another block diagram of the present 

invention for solving the Jc-dimensional assignment problem 
where 11 Jc" is greater than or equal to three (3) . 

Figure 8 is a flowchart describing the procedure for 
solving an n-dimensional assignment problem according to the 
2 0 second embodiment of the invention. 
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DETAILED DESCRIPTION OF THE IISTVENTION 
Referring now to the other figures in detail wherein like 
reference numerals indicate like elements throughout the 
several views, Figure 3 illustrates a system generally 
5 designated 100 for implementing the present invention. System 
100 comprises, for example, a radar station 102 (note sonar, 
microwave, infrared and other radiation bandwidths are also 
contemplated) for scanning a region which may be, for example, 
an aerial region (in aerial surveillance applications) or a 

10 work region (in robotic applications) and generating signals 
indicating locations of objects within the region. The 
signals are input to a converter 104 which converts the 
signals to data points or observations in which each object 
(or false positive) is represented by a single point. The 

15 output of the converter is input to and readable by a computer 
106. As described in more detail below, the computer 106 
assigns the points to respective tracks, and then displays the 
tracks and extrapolations of the tracks on a monitor 110. 
Also, the computer 106 determines an appropriate action to 

20 take based on the tracks and track extensions. For example, 
in a commercial application at an airport, the computer can 
determine if two aircraft being tracked are on a collision 
course and if so, signal a transmitter 112 to warn each 
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aircraft , or if a scheduled take-off will pose the risk of 
collision, delay the take-off. For a military application on 
a ship or base, the computer can determine subsequent 
coordinates of enemy aircraft and send the coordinates to an 
5 antiaircraft gun or missile 120 via a communication channel 
122 . In a robotic application, the computer controls the 
robot to work on the proper object or portion of the object. 

The invention generates .k-dimensional matrices where k is 
the number of images or sets of observation data in the look 

10 ahead window plus one. Then, the invention formulates a k- 
dimensional assignment problem as in (0.1) above. 

The ^-dimensional assignment problem is subsequently 
relaxed to a (k-1) -dimensional problem by incorporating one 
set of constraints into the objective function using a 

15 Lagrangian relaxation of this set. Given a solution of the 
(k-1) -dimensional problem, a feasible solution of the k- 
dimensional problem is then reconstructed. The (k-1) - 
dimensional problem is solved in a similar manner, and the 
process is repeated until it reaches the two-dimensional 

2 0 problem that can be solved exactly. The ideas behind the 
Lagrangian relaxation scheme are outlined next. 
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Consider the integer programming problem 



Minimize 



v(z) = c T z 



Subject to: 



Az < b, 
Bz < d, 

z ± is an integer for i £ I, 



(0.2) 



5 



where the partitioning of the constraints is natural in some 
sense. Given a multiplier vector u > 0, the Lagrangian 
10 relaxation of (0.2) relative to the constraints Bz < d is 
defined to be 



If the constraint set Bz < d is replaced by Bz = d, the 

non-negativity constraint on u is removed. S?=c T z+u T (Bz-d) is 

the Lagrangian relative to the constraints Bz < d, and hence 
2 0 the name Lagrangian relaxation. The following fact gives the 

relationship between the objective functions of the original 

and relaxed problems. 

FACT A.l. If z is an optimal solution to (0.2), then 

<P(u)< v(F) for all u > 0. If an optimal solution z of (0.3) 
25 is feasible for (0.3) , then z is an optimal solution for (0.2) 

and 0(u) - v(z) . 



<P(u) - Minimize 



{c T z+u T (Bz-d) } 



(0.3) 



Subject to: 



Az < b, 

z i is an integer for i G I. 
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This leads to the following algorithm: 
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Algorithm. Construct a sequence of multipliers {u k } k = 0 
converging to the solution u of Maximize {<P{u) : u > 0} and a 
corresponding sequence of feasible solutions {z k } k=0 of (0.2) as 
follows : 

5 1. Generate an initial approximation u 0 . 

2. Given u k , choose a search direction s k and a search 
distance a k so that <P{u k +a k s k ) > <P(u k ) . Update the multiplier 
estimate by u k+1 - u k + cx k s k . 

3. Given u k+1 and a feasible solution z k+1 (u k+1 ) of (0.3), 
10 recover a feasible solution z k+1 {u k+1 ) of the integer programming 

problem (0.2) . 

4. Check the termination criteria. If the algorithm is 
not finished, set k=k+l and return to Step 2. Otherwise, 
terminate the algorithm. 

15 If z is an optimal solution of (0.2), then 

0(u k ) < 0(u) < v(z) < v(z k ). 
Since the optimal solution £ = z(u) of (0.3) is usually not a 
feasible solution of (0.2) for any choice of the multipliers, 
@(u) is usually strictly less than vCz). The difference v(z)- 

20 0(u) is called the "duality gap," which can be estimated by 
comparing the best relaxed and recovered solutions computed in 
the course of maximizing 0(u). 
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Thus, to eliminate the 11 <" in the constraint formulations 
(e.g., (0.1) (a) - (0.1) (d) ) that resulted from false 
positives, the constraint sets have been modified, as 
indicated above, to incorporate a gap filler in each 
5 observation set to account for false positives. Thus, a zero 
for an index i p (1 < p < k) in i in problem (0.4) below 

indicates that the track extension represented by z*[ ± 
includes a false positive in the p th observation set. Note 
that this implies that a hypothesis be formed incorporating an 
10 observation with Jt-1 gap fillers, e.g., z 0 . . .o± p o.. .o* ^ P ^°- Thus, 
the resulting generalization of Problem Formulation (0.1) 
without the "less than" complication within the constraints is 
the following Jc-Dimensional Assignment Problems in which k >3 : 
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zJ . . . zJ c. . Z . 



VO i k =o 



Subject to: /J . . . Tj z. . = 1, i 2 = 1, . . -JV 2/ 



AT. , 



i —C\ 4 = A V = A -? = A 1 * " * .fc 



ij-0 i Jtl =0 i k =0 



for ij = 1, . . .Nj 
and j=2, . . .k - 1, 



z2 • ■ • 12 z i ± =1, i k = 1, . . .N k , 

Z i 1 ...i k E W n , 22 = 1,...,*, 



10 where c and z are similarly dimensioned matrices representing 
costs and hypothetical assignments. Note, in general, for 
tracking problems these matrices are sparse. 

After formulating the Jc-dimensional assignment problem as 
in (0.4), the present invention solves the resulting problem 

15 so as to generate the outputs required by devices 110, 112, 
122 and 130. For each observation set 0 ± received from 
converter 104 at time t ± where 1=1,...,^, the computer 
processes Oi in a batch together with the other observation 
sets Oi_ k+1/ . . . ,Oi and the track T^, i.e., Tj is the set of all 

2 0 tracks that have been defined up to but not including Oj . 
(Note, bold type designations refer to the vector of elements 
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for the indicated time, i.e., the set of all observations in 
the scan or tracks existing at the time, etc.) The result of 
this processing is the new set of tracks Ti_ k+1 and a set of 
cost weighted possible solutions indicating how the tracks 
5 might extend to the current time t ± . At time fc i+1 the batch 
process is repeated using the newest observation set and 
deleting the oldest. Thus, there is a moving window of 
observation sets which is shifted forward to always 
incorporate the most recent observation set. The effect is 

10 that input observation sets are reused for (k-1) cycles and 
then on the observation set's k th reuse each observation within 
the observation set is integrated into a track. 

Figure 7 illustrates various processes implemented upon 
receipt of each observation set. Except for the addition of 

15 the k- dimensional assignment solving process 30 0 and the 
modification to scoring process 154 to build data structures 
suitable for process 3 00, all processes in Figure 7 are based 
upon prior art. The following processes 150 and 152 extend 
previously defined tracks based on new observations. Gate 

2 0 formulation and output process 156 determines, for each of the 
previously defined tracks, a zone wherein the track may 
potentially extend based on limits of velocity, 
maneuverability and radar precision. One such technique to 



28 



CSUR.01USRI 

accomplish this is the cone method described previously. The 
definition of the zone is passed to gating process 150. When 
a new observation set O t is received, the gating process 150 
will match each member observation with the zone for each 
5 member of the hypothetical set . After all input 

observations from O i are processed, the new hypothesis set hi 
is generated by extending each track of the prior set of 
hypothetical tracks either with missed detect gap fillers 
or with all new observation elements satisfying the track's 

10 zone. This is a many-to-many matching in that each hypothesis 
member can be extended to many new observations and each new 
observation can be used to extend many hypotheses. It, 
however, is not a full matching in that any hypothesis will 
neither be matched to all observations nor vice versa. It is 

15 this matching characteristic that leads to the sparse matrices 
involved in the tracking process. Subsequently, gating 150 
forwards the new hypothesis set h ± to filtering process 152 . 
Filtering process 152 determines a smooth curve for each 
member of h 2 . Such a smooth curve is more likely than a sharp 

2 0 turn from each point straight to the next point. Further, the 
filtering process 152 removes small errors that may occur in 
generating observations. Note that in performing these tasks, 
the filtering process 152 preferably utilizes a minimization 
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of a least squares test of the points in a track hypothesis or 
a Kalman filtering approach. 

As noted above, the foregoing track extension process 
requires knowledge of a previous track. For the initial 
5 observations, the following gating process 158 and filtering 
process 160 determine the "previous track" based on initial 
observations. In determining the initial tracks, the points 
from the first observation set form the beginning points of 
all possible tracks . After observation data from the next 

10 observation set is received, sets of simple two-point straight 
line tracks are defined. Then, promotion, gate formulation, 
and output step 162 determines a zone in which future 
extensions are possible. Note that filtering step 160 uses 
curve fitting techniques to smooth the track extensions 

15 depending upon the number of prior observations that have 
accumulated in each hypothesis. Further note that promotion, 
gate formulation and output process 162 also determines when 
sufficient observations have accumulated to form a basis for 
promoting the track to processes 150 and 152 as described 

2 0 above . 

The output of each of the filtering processes 152 and 160 
is a set of hypothetical track extensions. Each such 
extension contains the hypothetical initial conditions (from 
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the previous track) , the list of observations incorporated in 
the extension, and distance between each observation and the 
filtered track curve. Scoring process 154 determines the 
figure of merit or cost of an observation being an extension 
5 of a track. In one embodiment, the cost is based on the 
above-mentioned distance although the particular formula for 
determining the cost is not critical to the present invention. 
A preferred formula for determining the cost utilizes a 
negative log likelihood function in which the cost is the 

10 negative of the sum of: (a) the logs of the distances 
normalized by sensor standard deviation parameters, and (b) 
the log likelihoods for events related to track initiation, 
track termination, track maneuver, false negatives and false 
positives. Note that track maneuvers are detected by 

15 comparing the previous track curve with the current extension. 
Further note that some of the other events related to, for 
example, false negatives and false positives are detected by 
analyzing the relative relationship of gap fillers in the 
hypothesis. Thus, after determining that one of these events 

2 0 occurred, a cost for it can be determined based upon suitable 
statistics tables and system input parameters. The negative 
log likelihood function is desirable because it permits 
effective integration of the useful components. Copies of the 
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set of hypothetical track extensions which are scored are 
subsequently passed directly to one of the gate formulation 
and output steps 156 and 162. Note that the scoring process 
154 also arranges the actual scores in a sparse matrix based 
5 upon observation identifiers , and passes them to Jc-dimensional 
assignment problem solving process 300. 

The assignment solving process 3 00 is described below. 
Its output is simply the list of assignments which constitute 
the most likely solution of the problem described by Equation 

10 (0.4). Note that both gate formulation and output processes 
156 and 162 use (at different times) the list of assignments 
to generate the updated track history T ± to eliminate or prune 
alternative previous hypotheses that are prohibited by the 
actual assignments in the list, and, subsequently, to output 

15 any required data. Also note that when one of the gate 
formulation and output processes 156 and 162 accomplish these 
tasks, the process will subsequently generate and forward the 
new set of gates for each remaining hypothesis and the 
processes will then be prepared to receive the next set of 

2 0 observations. In one embodiment, the loop described here will 
generate zones for a delayed set of observations rather than 
the subsequent set. This permits processes 156 and 162 to 
operate on even observation sets while the scoring step 154 
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and Tc-dimerisional solving process operate on odd sets of 
observations, or vice versa. 

The assignment solving process 300 permits the present 
invention to operate with a window size of dimension k-1 for 
5 some k > 3 . The upper limit on k depends only upon the 
computational power of the computer 106 and the response time 
constraints of system 100. The k-1 observation sets within 
the processing window plus the prior track history result in 
a Jt-dimensional Assignment Problem as described by Problem 
10 Formulation (0.4). The present invention solves this 
generalized problem including the processes required to 
consider false positives and negatives, and also the processes 
required to consider sparse matrix problem formulations. 

15 I . A FIRST EMBODIMENT OF THE 

^-DIMENSIONAL ASSIGNMENT SOLVER 300 

In describing a first embodiment of the k- dimensional 

assignment solver 3 00, it is worthwhile to also discuss the 

process of Figure 4 which is used by the solver 300. Figure 

2 0 4 illustrates use of the Frieze and Yadagar process as shown 

in prior art for transforming a three-dimensional assignment 

problem into a two-dimensional assignment problem and then use 

a hill climbing algorithm to solve the three-dimensional 

assignment problem. The solver 3 00 uses a Lagrangian 
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Relaxation technique (well known in the art) to reduce the 
dimension of an original ^-dimensional assignment problem 
(k > 3) down to a three-dimensional problem and then use the 
process of Figure 4 to solve the three-dimensional problem. 
5 Further note that the Lagrangian Relaxation technique is also 
utilized by the process of Figure 4 and that in using this 
technique the requirement that each point is assigned to one 
and only one track is relaxed. Instead, an additional cost, 
which is equal to a respective Lagrangian Coefficient u, is 

10 added to the cost or objective function (0.1) (a) whenever a 
point is assigned to more than one track. This additional 
cost can be picked to weight the significance of each 
constraint violation differently, so this additional cost is 
represented as a vector of coefficients u which are correlated 

15 with respective observation points. Hill climbing will then 
develop a sequence of Lagrangian Coefficients sets designated 
(u 0 , . . . ,u jf u j+1/ . . . ,Up) . That correspond to an optimum solution 
of the two-dimensional assignment problem. The assignments at 
this optimum solution are then used to "recover" the 

2 0 assignment solution of the three-dimensional assignment 
problem. 

In step 200 of Figure 4, initial values are selected for 
the u 0 coefficients. Because the Lagrangian Relaxation process 
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is iterative, the initial values are not critical and are all 
initially selected as zero. In step 202, these additional 
costs are applied to the objective function (0.1) (a). With 
the addition of the costs u, the goal is still to assign the 
5 points which minimize the total cost. This transforms 
Equation (0.1) (a), written for k= 3 and altered to exclude 
mechanisms related to false positives and negatives, into 
objective function (1.1) (a) . In the first iteration it is not 
necessary to consider the u matrix because all u values are 

10 set to zero. To relax the requirement that each point be 
assigned to one and only one track, the constraint Equation 
(0.1) (d) is deleted, thereby permitting points from the last 
image to be assigned to more than one track. Note that while 
any axis can be chosen for relaxation, observation constraints 

15 are preferably relaxed. The effect of this relaxation is to 
create a new problem which must have the same solution in the 
first two axes but which can have a differ solution in the 
third axis. The result is constraints (1.1) (b-d) . 
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W, N, 



(a) Minimize }1 E ^( c i ± ± -u.))z 

i 1= o i 2 =o 1 2 3 3 1 2 3 

AT, N, 



(b) Subject to: 



i 2 =l i 3 =l 1 2 3 



(c) £ J) z. , . <1, i 2 = 1,...,^, 

i 3 =l 1 2 3 

(d) ^1,1,1* <0,1> v 



z . 



Step 2 04 then generates from the three-dimensional problem 
described by Problem Formulation (0.1) a new two-dimensional 
problem formulation which will have the same solution for the 

10 first two indices. As Problem Formulation (1.1) has no 
constraints on the 3 rd axis, any value within a particular 3 rd 
axis can be used in a solution, but using anything other than 
the minimum value from any 3 rd axis has the effect of 
increasing solution cost. Conceptually, the effect of step 

15 2 04 is to change the three-dimensional arrays in Problem 
Formulation (1.1) into two-dimensional arrays as shown in 
Problem Formulation (1-2) and to generate the new two- 
dimensional matrix m ±li2 defined as shown in Equation (1.3) . 
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[a) Minimize 



(b) Subject to: 



(c) 
(d) 



i^O i 2 =0 



Min: 



(c. . . -u. ) 



V i, 



i 2 -i i 3 =i ^ 



M, N, 



j£ z^^l, i 2 = 1,...,N 2 , 

i,=i i 3 =i 1 3 

z. e {0,1} V z. , 



(1.2) 



m. . = Minimum { c . . - u . I i = l, . . . , AT }■ 



(1.3) 



The cost or objective function for the reduced problem as 
defined by (1.2) (a), if evaluated at all possible values of u 
is a surface over the domain of u. This surface is referred 

10 to as 0(u) and is non-smooth but provable convex (i.e., it has 
a single peak and several other critical characteristics which 
form terraces) . Due to the convex characteristics of <P(u) , 
the results from solving Problem Formulation (1.2) at any 
particular Uj can be used to generate a new set of coefficients 

15 u j+1 whose corresponding cost value is closer to the peak of 
@(u) . The particular set of Lagrangian Coefficients that will 
generate the two-dimensional problem resulting in the maximum 
cost is designated u p . To recover the solution to the three- 



37 



CSUR.01USRI 

dimensional assignment problem requires solving the Equation 
(1.2) problem corresponding to u p . 

In step 2 06, the two-dimensional problem is solved 
directly using a technique known to those skilled in the art 
5 such as Reverse Auction for the corresponding cost and 
solution values. This is the evaluation of one point on the 
surface or for the first iteration &(u 0 ) . 

Thus, after this first iteration, the points have been 
assigned based on all "u" values being arbitrarily set to 

10 zero. Because the "u" values have been arbitrarily assigned, 
it is unlikely that these assignments are correct and it is 
likely that further iterations are required to properly assign 
the points. Step 208 determines whether the points have been 
properly assigned after the first iteration by determining if 

15 for this set of assignments whether a different set of "u" 
values could result in a higher total cost. Thus, step 2 08 is 
implemented by determining the gradient of objective function 
(1.2) (a) with respect to Uj . If the gradient is substantially 
non-zero (greater than a predetermined limit) then the 

2 0 assignments are not at or near enough to the peak of the 
surface @(u) (decision 210) , and the new set of Lagrangian 
Coefficients u j+1 is determined. 
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Hill climbing Step 212 determines the u j+1 values based 
upon the Uj values, the direction resulting from protecting the 
previous gradient into the U 3 domain, and a step size. The 
solution value of the two-dimensional problem is the set of 
5 coefficients that minimize the two-dimensional problem and the 
actual cost at the minimum. Those coefficients augmented by 
the coefficients stored in m il±2 permit the evaluation (but not 
the minimization) of the cost term in Problem Formulation 

(1.1) . These two cost terms are lower and upper bounds on the 
10 actual minimized cost of the three-dimensional problem, and 

the difference between them in combination with the gradient 
is used to compute the step size. 

With this new set of u values, steps 202-210 are repeated 
as a second iteration. Steps 212 and 202-210 are repeated 
15 until the gradient as a function of u determined in step 2 08 
is less than the predetermined limit. This indicates that the 
u p values which locate the peak area of the <P{u) surface are 
determined and that the corresponding Problem Formulation 

(1.2) has been solved. Step 214 will attempt to use the 
2 0 assignments that resulted from this particular two-dimensional 

assignment problem to recover the solution of the three- 
dimensional assignment problem as described below. If the 
limit was chosen properly so that the u values are close 
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enough to the peak, this recovery will yield the proper set of 
assignments that rigidly satisfies the constraint that each 
point be assigned to one and only one track. However, if the 
u values are not close enough to the peak, then the limit 
5 value for decision 210 is reduced and the repetition of steps 
212 and 202-210 is continued. 

Step 214 recovers the three-dimensional assignment 
solution by using the assignment values determined on the last 
iteration through step 208. Consider the two-dimensional z 

10 assignment matrix to have l ! s in the locations specified by 
the list L 1 ={a i , b ± ) N i=lm If the three-dimensional z matrix is 
specified by placing 1 1 s at the location indicated by the list 
L 2 ={a i/ b if ta ah )^ =1 then the result is a solution of Problem 
Formulation (1.1). Let L 3 = {m aib .) f ml be the list formed by the 

15 third index. If each member of L 3 is unique then the L 2 
solution satisfies the third constraint so it is a solution to 
Problem Formulation (0.1). When this is not the case, 
recovery determines the minimal substitutions required within 
list L 3 so that it plus L x will be a feasible solution, i.e., 

20 a solution which satisfies the constraints of a problem 
formulation, but which may not optimize the objective function 
of the problem formulation. This stage of the recovery 
process is formulated as a two-dimensional assignment problem 
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as follows. Form a new cost matrix [c ifj ]*[ /j9l where c lfj =c a±bifJ for 
j=l...N ± and the N ± term is the total number of cost elements 
in the selected row of the three-dimensional cost matrix. 
Attempt to solve this two-dimensional problem for the minimum 
5 using two constraints sets. If a feasible solution is found 
then the result will have the same form as list L x . Replace 
the first set of indices by the indicated {a if b ± ) pairs taken 
from list L 2 and the result will be a feasible solution of 
Problem Formulation (0.4) . If no feasible solution to the new 
10 two-dimensional problem exists then further effort to locate 
the peak of cP(u) is required. 



I . 1 . Generalization to a Multi-Dimensional 
Assignment Solving Process 

15 Let M be a fixed integer and assume that M is the 

dimension of the initial assignment problem to be solved. 
Thus, initially, the result of the scoring step 154 is a M- 
dimensional Cost Matrix which is structured as a sparse matrix 
(i.e., only a small percentage of the entries in the cost and 

2 0 assignment matrices are filled or non-zero) . Individual cost 
elements represent the likelihood that a track T ± as extended 
by the set of observations {O^j i = l,...,M-l}, is not valid. 
Because the matrix is sparse the list of cost elements is 
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stored as a packed list, and then for each dimension of the 
matrix, a vector of a variable length list of pointers to the 
cost elements is generated and stored. This organization 
means that for a particular observation 0 ±j for the j th list in 
5 the 2 th vector will be a list of pointers to all hypotheses in 
which 0 ±j participates. This structure is further explained in 
the following section dealing with problem partitioning. 

The objective of the assignment solving process is to 
select from the set of all possible combinations of track 

10 extensions a subset that satisfies two criteria. First, each 
point in the subset of combinations should be assigned to one 
and only one track and therefore, included in one and only one 
combination of the subset, and second, the total of the 
scoring sums for the combinations of the subset should be 

15 minimized. This yields the following M-dimensional equations 
where k = M: 
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(a) Minimize v k (z k )= 22 ■ • • 22 c i ...i z i ...i (1.4) 

(b) Subject to: • • • z i i = 1* i 2 = 1, . . .N lf 

i 2 =o vo 



(c) 



w. 1 n.^ 



5^ ••• ^ ^ • • • 5^ ^...i. 1' 

i^O v^o i j+1 = 0 1^0 

for i . = 1 , . • . , N . and j = 2 , . . . , k-1, 



10 (d) X, • • • X, ^...i^ 1 ' i j t = 1 / ' ' - 



(e) z^^e {0, 1} for all i lf . . . f l k , 



and where c k is the cost matrix [c^.^^] which is a function of 
15 the distance between the observed point z k and the smoothed 
track determined by the filtering step, and v k is the cost 
function. This set of equations is similar to the set 
presented in Problem Formulation (0.4) except that it includes 
the subscript and superscript k notation. Thus, in solving 
2 0 the M-dimensional Assignment Problem the invention reduces 
this problem to an (M-l) -dimensional Assignment Problem and 
then to an (N-2 ) -dimensional Assignment Problem, etc. 
Further, the symbol ke{3, . . . ,M} customizes Problem Formulation 
(1.4) to a particular relaxation level. That is, the notation 
25 is used to reference data from levels relatively removed as in 
c k+1 are the cost coefficients which existed prior to this level 
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of relaxed coefficients c k . Note that actual observations are 
numbered from 1 to W i; where N ± is the number of observations 
in observation set i. Further note that the added 0 
observation in each set of observations is the unconstrained 
5 "gap filler." This element serves as a filler in substituting 
for missed detects, and a sequence of these elements including 
only one true observation represents the possibility that the 
observation is a false positive. Also note that by being 
unconstrained a gap filler may be used in as many hypotheses 

10 as required. 

While direct solution to (1.4) would give the precise 
assignment, the solution of Jc-dimensional equations directly 
for large k is too complex and time consuming for practice. 
Thus, the present invention solves this problem indirectly. 

15 The following is a short description of many aspects of 

the present invention and includes some steps according to the 
prior art. The first step in solving the problem indirectly 
is to reduce the complexity of the problem by the previously 
known and discussed Lagrangian Relaxation technique. 

2 0 According to the Lagrangian Relaxation technique, the absolute 
requirement that each point is assigned to one and only one 
track is relaxed such that for some one image, points can be 
assigned to more than one track. However, a penalty based on 
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a respective Lagrangian Coefficient u k is added to the cost 
function when a point in the image is assigned to more than 
one track. The Lagrangian Relaxation technique reduces the 
complexity or "dimension" of the formulation of the assignment 
5 problem because constraints on one observation set are 
relaxed. Thus, the Lagrangian Relaxation is performed 
iteratively to repeatedly reduce the dimension until a two- 
dimensional penalized cost function problem results as in 
Problem Formulation (1.1). This two-dimensional problem is 

10 solved then directly by a previously known technique such as 
Reverse Auction. The penalized cost function for the two- 
dimensional problem defines a valley or convex shaped surface 
which is a function of various sets of {u k \k=3 , . . . , M} penalty 
values and one set of assignments for the points in two 

15 dimensions. That is, for each particular u 3 there is a 
corresponding two-dimensional penalized cost function problem 
and its solution. Note that the solution of the two- 
dimensional penalized cost function problem identifies the set 
of assignments for the particular u 3 values that minimize the 

20 penalized cost function. However, these assignments are not 
likely to be optimum for any higher dimensional problem 
because they were based on an initial arbitrary set of u k 
values. Therefore, the next step is to determine the optimum 
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assignments for the related three-dimensional penalized cost 
function problem. There exists a two-dimensional hill shaped 
function 0 which is a graph of the minimums of all penalized 
cost functions at various sets of assignments in two 
5 dimensions. For the three-dimensional problem, the function 
0 can be defined based on the solution to the foregoing two- 
dimensional penalized cost function. By using the current u 3 
values and the {u k |Jc > 3} values originally assigned, the 
gradient of the hill -shaped function <P is determined, which 

10 points toward the peak of the hill. By using the gradient and 
u 3 values previously selected for the one point on the hill 
(corresponding to the minimum of the penalized cost function 
0 ) as a starting point, the u 3 values can be found for which 
the corresponding problem will result in the peak of the 

15 function 0. The solution of the corresponding two-dimensional 
problem is the proper values for two of the three sets of 
indices in the three-dimensional problem. These solution 
indices can select a subsection of the cost array which maps 
to a two-dimensional array. The set of indices which minimize 

20 the two-dimensional assignment problem based on that array 
corresponds to the proper assignment of points in the third 
dimension. The foregoing "recovery" process was known in the 
prior art, but it is modified here to adjust for the sparse 
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matrix characteristic. The next task is to recover the 
solution of the proper u 4 values for the four-dimensional 
problem. The foregoing hill climbing process will not work 
again because the foregoing hill climbing process when 
5 required to locate the four-dimensional u 4 values for the peak 
of &* requires the exact definition of the function & (as was 
available in the case of $?) or an always less than 
approximation of &* , whereas the iteration can result in a 
greater than approximation of . According to the present 

10 invention, another three-dimensional function ¥ is defined 
which is a "less than approximation" the three-dimensional 
function <P and which can be defined based on the solution to 
the two-dimensional penalized cost function and the previously 
assigned and determined u k values. Next, the gradient of the 

15 function ¥ is found and hill climbing technique used to 
determine the u 4 values at the peak. Each selected u 4 set 
results in a new three-dimensional problem and requires the 
two-dimensional hill climbing based upon new two-dimensional 
problems. At the peak of the three-dimensional function ¥, 

20 the solution values are a subset of the values required for 
the four-dimensional solution. Recovery processing extends 
the subset to a complete solution. This process is repeated 
iteratively until the u k values that result in a corresponding 
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solution at the peak of the highest order functions ¥ and 0 
are found. The final recovery process then results in the 
solution of Jc-dimensional problem. 

Figure 5 illustrates process 300 in more detail. 

5 

1.2. Problem Formulation 
In problem formulation step 310, all data structures for 
the subsequent steps are allocated and linked by pointers as 
required for execution efficiency. The incoming problem is 

10 partitioned as described in the subsequent section. This 
partitioning has the effect of dividing the incoming problem 
into a set of independent problems and thus reducing the total 
workload. The partitioning process depends only on the actual 
cost matrix so the partitioning can and is performed for all 

15 levels of the relaxation process. 

1.2.1, Relaxation and Recovery 
Step 32 0 begins the Lagrangian Relaxation process for 
reducing the M-dimensional problem by selecting all Lagrangian 
2 0 Coefficient u M penalty values initially equal to zero. The 
Lagrangian Coefficients associated with the constraint set 
are an (2^+1) -element vector. The reduction of this M- 
dimensional problem in step 322 to a (Af-1) -dimensional problem 
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uses the two step process described above. First, a penalty 
based on the value of the respective u M coefficient is added to 
the cost function when a point is assigned to more than one 
track and then the resultant cost function is minimized. 
5 However, during the first iteration, the penalty is zero 
because all u M values are initially set to zero. Second, the 
requirement that no point from any image can be assigned to 
more than one track is relaxed for one of the images. In the 
extreme this would allow a point from the relaxed image to be 

10 associate with every track. However, the effect of the 
previous penalty would probably mean that such an association 
would not minimize the cost. The effect of the two steps in 
combination is to remove a hard constraint while adding the 
penalty to the cost function so that it operates like a soft 

15 constraint. For step 322 this two-step process results in the 
following penalized cost function problem with k-M\ 
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(a) <P k (u k ) = Minimize: (p k (u k f z k ) 



Z^t • • • Z^ . .1. z i,. . .2. Z-f U 2\ 

i 1= 0 i k =0 i k =0 

1^0 2^=0 1^=0 



V 0 ijt-i 



(b) Subject to: 



(c) 



(1.5) 



(d) 



Z^ * * * Z^f Z*/ • " • Z^ ^1-...!^. i 

i x = 0 i j+1 = 0 i k =o 

for ij = l f . . .Nj 
and j-2, . . .k-1, 

z i 1 ...i k e fa* 1 } for all i lm ..l k . 



10 Because the constraint on single assignment of elements from 
the last image has been eliminated, an (Af-1) -dimensional 
problem can be developed by eliminating some of the possible 
assignments. As shown in Equations (1.6), this is done by 
selecting the smallest cost element from each of the M th axis 

15 vectors of the cost matrix. Reduction in this manner yields 
a new, lower-order penalized cost function defined by 
Equations (1.6) which has the same minimum cost as does the 
objective function defined by (1.5) (a) above. 

The cost vectors are selected as follows. Define a new 



2 0 cost array c^ mmi by: 
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c^^^ Minimum {c^.. ifc - 1^= 0, 1, . . . , Z^j, (1.6) 
for {i lf . . . , i k ^) * (0, . . . , 0) , 

The resulting (Af-1) -dimensional problem is (where k=M) : 
tfkfu*; = Minimize: v k _ 1 (z k ' 1 )= X) • • • X, c i\?. .i^ ^i,?. .1, 2 

Subject to: X, • • • $J z i~ 1 i = 1 ' ^ = 2 ' * • - N i' ^ ■ 7 ^ 

i 2 =o vo 



X/ • • • S ]L/ • • • iLf z i 1 . . .i k 2 ' 
i 1= 0 i J+1 = 0 i k =0 

for ij = 1, . . .Nj 
and j=2, . . .k-1, 

2-/ . . .1. ~ 1 ' -^Jc-l ™ i/''*^k-l/ 

e {0,1} for all i,...!*.,. 



10 Assignment Problem (1.4) and Problem Formulation (1.7) differ 
only in the dimension M vs. M-l, respectively. An optimal 
solution to (1.7) is also an optimal solution to equation 
(1.5). This relationship is the basis for an iterative 
sequence of reductions indicated by steps 320-332 through 330- 

15 332 and 200-204 in which the penalized cost function problem 
is reduced to a two-dimensional problem. As these formula 
will be used to describe the processing at all levels, the 
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lowercase k is used except where specific reference to the top 
level is needed. In step 206, the two-dimensional penalized 
cost function is solved directly by the prior art Reverse 
Auction technique. Each execution of 206 produces two 
5 results, the set of z 2 values that minimize the problem and the 
cost that results v 2 when these z values are substituted into 
the objective function (1.7) (a) . 

In step 2 08, according to the prior art, solution z 2 
values are substituted into the two-dimensional derivative of 

10 the surface <t> 2 . The result indicates how the value of u 3 
should be adjusted so as to perform the hill climbing 
function. As was previously described the objective is to 
produce a sequence of u? values which ends when the value is 
in the domain of the peak of the surface 0. The section 

15 "Determining Effective Gradient" describes how new values are 
computed and how it is determined that the current points to 
the peak of @ 2 „ When no further adjustment is required the 
flow moves to step 214 which will attempt to recover the 
three-dimensional solution as previously described. When 

2 0 further adjustment is required then the flow progresses to 
step 212 and the new values of u k are computed. At the two- 
dimensional level the method of the prior art could be used 
for the hill climbing procedure. However, it is not practical 
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to use this prior art hill climbing technique to determine the 
updated Lagrangian Coefficients u k or the Max on the next (or 
any) higher order surface <P because the next (or any) higher 
dimensional function 0 cannot be defined based on known 
5 information. 

Instead, the present invention defines a new function 
based on known information which is useful for hill climbing 
from the third to the fourth and above dimensions, i.e., u k 
values which result in z values that are closer to the proper 

10 z values for the highest k-dimension. This hill climbing 
process (which is different than that used in the prior art of 
Figure 4 for recovering only the three-dimensional solution) 
is used iteratively at all lower dimensions k of the M- 
dimensional problem (including the three-dimensional level 

15 where it replaces prior art) even when k is much larger than 
three. Figure 6 helps to explain this new hill climbing 
technique and illustrates the original k-dimensional cost 
function v k of Problem Formulation (1.4) . However, the actual 
k-dimensional cost surface v k defined by (1.4) (a) comprises 

2 0 scalar values at each point described by k-dimensional vectors 
and as such can not be drawn. Nevertheless, for illustration 
purposes only, Figure 6 ignores the reality of ordering 
vectors and illustrates a concave function v fo (z k ) to represent 
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Equation (1.4) . The surface v k is illustrated as being smooth 
to simplify the explanation although actually it can be 
imagined to be terraced. The goal of the assignment problem 
is to find the values of "z k ; these values minimize the k- 
5 dimensional cost function v k . 

For purposes of explanation, assume that in Figure 6, Jc=4 
(the procedure is used for all k >3) . This problem is reduced 
by two iterations of Lagrangian Relaxation to a two- 

(k-l) i 

dimensional penalized cost function 0^ . This cost 
10 function, and all other cost functions described below, are 
also non- smooth and continuous but are illustrated in Figure 

(k-l) x 

6 as smooth for explanation purposes. Solving the 0. 
problem results in one set of z 2 assignments and the value of 

(k~l) i (k-l) 1 

@(k-i)± at the point . A series of functions <pj , . . . , <p 1 
15 each generated from a different u 3 are shown. The particular 
series illustrates the process of locating the peak of i > {k - 1)i ^ 

(k-l) i (k-l) i 

The two-dimensional penalized cost functions <pj , . . . , <p x 
can be solved directly. Each such solution provides the 
information required to calculate the next u 3 value. Each 
20 iteration of the hill climbing improves the selection of u 3 
values, i.e., yields 0 2 problem whose solution is closer to 
those at the solution of the 0 3 problem. The result of solving 

(k-l). 

0 1 is values that are on both < P (k ^ 1)i and & k . Figure 6 
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illustrates the surface & k which comprises the minimums of all 
^-dimensional penalized cost function surfaces <p kf i.e., if the 
Problem Formulations (1.5) and (1.7) were solved at all 
possible values of u k the function @ k (u k ) would result. The 
5 surface 0 k is always less than the surface v*. except at the 
peak as described in Equation (1.8) and its maximum occurs 
where the surface v k is minimum. Because the surface & k 
represents the minimum of the surface <P k , any point on the 
surface & k can mapped to the surface v k . The function 0 k 
10 provides a lower bound on the minimization problem described 
by (1.4) (a) . Let ~z k be the unknown solution to Problem 
Formulation (1.4) and note that: 

® k [u k )< v k {z k )z v k _ x {z k ) . (1.8) 

15 Consequently, the values z k at the peak of 0 k (i.e., the 
greatest lower bound on the cost of the relaxed problem) , can 
be substituted into the Jc-dimensional penalized cost function 
to determine the proper assignments. Consequently, the 
present invention attempts to find the maximum on the surface 

20 0 k . However, it is not possible to use the prior art hill 
climbing to hill climb to the peak of & k because the definition 
of 0 k requires exact knowledge of lower order functions &. As 
the solution of <p\ is not the exact solution, in that higher 
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order values of u are not yet optimized, its value can be 
larger than the true peak of & k . As such it is not a lower 
bound on v k and it can not be used to recover the required 
solution . 

5 Instead, the present invention defines all auxiliary 

functions W k which are based only on the solution to the two- 
dimensional penalized cost function problem, lower order 
values z k and values u k determined previously by the reduction 
process. The function W k is a less than approximation of <P kr 

10 and its gradient is used for hill climbing to its peak. The 
values z k at the peak of the function W k are then substituted 
into Problem Formulation (1-5) to determine the proper 
assignments. To define the function W kf the present invention 
explicitly makes the function <P k (u k ) a function of all higher 

15 order sets of Lagrangian Coefficients with the expanded 
notation: <P k (u k ; u k+1 , . . . ,u k ) . Then, a new set of functions W, 
is defined recursively, using the ^'s domain: 

Vu 3 ) = ^ 2 + X, =$ 3 (u 3 ;u 4 ,...,u K ), (1.9) 

20 

where v 2 is the solution value for the most recent two- 
dimensional penalized cost function problem. For k >3 
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WAu 3 , . . .,u* _1 ;u*) = 



®Au k ;u k+1 ,...,u K ), 



if known. 



^ 1 (u 3 f ... f u k ' 2 ;u k - 1 ) + }2 (1.10) 



otherwise. 



From the definitions of <P k and (Problem Formulation (1.7) 
compared with Problem Formulation (1.5)) : 

<b k ( u k ; u * +1 , . . . , u M ) = Vl ( z ^ ) + Jj 

it follows that : 

V" 3 ) = v 2+22 « 3 (u 3 ;uS . . .,u M )* v 3 (z 3 ) 
i 3 =o 

and with that Equation (1.8) is extended to: 
Vu 3 , • •.,u k - 1 ?u k )<® k (u k ?u k+1 , . ..,u M )<v k (u k )<v k (u k ) (1.11) 

5 This relationship means that either & k or W k may be used in 
hill climbing to update the Lagrangian Coefficients u. ^ is 
the preferred choice, however it is only available when the 
solution to Problem Formulation (1.5) is a feasible solution 
to Problem Formulation (1.4) (as in hill climbing from the 
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second to third dimension which is why prior art worked) . For 
simplicity in implementation, the function W k is defined so 
that it equals $ k when either function could be used. It is 
therefore always used, even for hill climbing from the second 
5 dimension. The use of the function W which is based on 
previously assigned or determined u values and not higher 
order u values which are not yet determined, is an important 
feature of the present invention. 

10 1.2.2. Determining Effective Gradient 

After the function W k is defined, the next steps of hill 
climbing/peak detection are to determine the gradient of the 
function ¥ kI determine an increasing portion of the gradient 
and then move up the surface W k in the direction of this 

15 increasing portion of the gradient. As shown in Figure 6, any 
upward step on the surface W kr for example, to the minimum of 
the 0 k j will yield a new set of values u k (to the "left") that 
is closer to the ideal set of values u k which correspond to the 
minimum of the function v k . While it is possible to 

20 iteratively step up this surface W k with steps of fixed size 
and then determine if the peak has been reached, the present 
invention optimizes this process by determining the single 
step size from the starting point at the minimum of 0 k . that 
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will jump to the peak and then calculate the values u k at the 
peak. Once the values u at the peak of ¥ k are determined, then 
the values u k can be substituted into Problem Formulation (1.5) 
to determine the proper assignment. (However, in a more 
5 realistic example, where k is much greater than three, then 
the values u k at the peak of the function ¥ along with the 
lower-order values u k and those assigned and yielded by the 
reduction steps are used to define the next higher level 

Q 

l?% function ¥. This definition of a higher order function ¥ and 

rp| 10 hill climbing process are repeated iteratively until ¥ kr the 
4;! peak of ¥ k , and the values u k at the peak of ¥ k are identified.) 

The following is a description of how to determine the 
Q gradient of each surface ¥ and how to determine the single 

TM step size to jump to the peak from the starting point on each 

W 15 surface ¥. 

As noted above, each surface ¥ is non-smooth. Therefore, 
if a gradient was taken at a single point, the gradient may 
not point toward the peak. Therefore, several gradients (a 
"bundle") are determined at several points on the surface ¥ k 
20 in the region of the starting point (i.e., minimum of <p kj ) and 
then averaged. Statistically, the averaged result should 
point toward the peak of the surface ¥ k . Wolfe's Conjugate 
Subgradient Algorithm (P. Wolfe. A method of conjugate 
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subgradients for minimizing non-dif f erentiable functions. 
Mathematical Programming Study, 3:145-173, 1975; P. Wolfe. 
Finding the nearest point in a polytope. Mathematical 
Programming Study, 11:12 8-14 9, 1976) for minimization was 
5 previously known in another environment to determine a 
gradient of a non- smooth surface using multiple subgradients 
and can be used with modification in the present invention. 
The modification to Wolfe's algorithm uses the information 
generated for ¥ k (u k3 , . . . ,u kk ~ 2 ; u kk_1 ) as the basis for calculating 
10 the subgradients. The definition of a subgradient v of 
^(u* 3 , . . . ,u kk ) is any member of the subdif f erential set defined 
as : 



^t(u)={v G R N * +1 | (?^(u k3 , . . . / u kk " 1 ;u / ) -^(u* 3 , . . . / u 3t " 1 ;u k ) ) 
15 >v T (u'-u k ) V u'E R N * +1 }, 

where v T is the transpose of v. 

Next, a subgradient vector is determined from this 
function. If z k is the solution of Problem Formulation (1.5), 
then differentiating ^(u 3 ,...,^" 1 ; u k ) with respect to u k k and 
2 0 evaluating the result with respect to the current selection 
matrix z k yields a subgradient vector: 
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g= (0, (i 



j!Li • • • 111 

ii=0 i^=0 



V 1 ' 



The iterative nature of the solution process at each dimension 
yields a set of such subgradients . Except for a situation 
described below, where the resultant averaged gradient does 
not point toward the peak, the most recent set of such 
5 subgradients are saved and used as the "bundle" for the peak 
finding process for this dimension. For example, at the 
level k there is a bundle of subgradients of the surface W k 
near the minimum of the surface <p k . determined as a result of 
solving Problem Formulation (1.4) at all lower levels. This 

10 bundle can be averaged to approximate the gradient . 
Alternately, the previous bundle can be discarded so as to use 
the new value to initiate a new bundle. This choice provides 
a way to adjust the process to differing classes of problems, 
i.e., when data is being derived from two sensors and the 

15 relaxation proceeds from data derived from one sensor to the 
other then the prior relaxation data for the first sensor 
could be detrimental to performance on the second sensors 
data . 
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1.2.3. Step Size and Termination Criterion 
After the average gradient of the surface W k is 
determined, the next step is to determine a single step that 
will jump to the peak of the surface ¥ k . The basic strategy 
5 is first to specify an arbitrary step size, and then calculate 
the value of ¥ k at this step size in the direction of the 
gradient. If the value of W k is larger than the previous one, 
this probably means that the step has not yet reached the 
peak. Consequently, the step size is doubled and a new value 

10 of W k is determined and compared to the previous one. This 
process is repeated until the new value of W k is less than the 
previous one. At that time, the step has gone too far and 
crossed over the peak. Consequently, the last doubling is 
rolled-back and that step size is used for the iteration. If 

15 the initial estimated ¥ k is less than previous value then the 
step size is decreased by 50%. If this still results in a 
smaller value of W kr then the last step is rolled back and the 
previous step size is decreased by 25%. The following is a 
more detailed description of this process. 

2 0 With a suitable bundle of subgradients determined as just 

described, Wolfe's algorithm can be used to determine the 
effective subgradient d and the upgraded value Uj +1 . From the 
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previous iteration, or from an initial condition, there exists 
a step length value t. The value, 

u + = + td 

is calculated as an estimate of u k +1 . To determine if the 
5 current step size is valid we evaluate ?k(u 3 , . . . ,u k ~ 2 ; u + ) . If 
the result represents an improvement then we double the step 
size. Otherwise we halve the step size. In either case a new 
u + is calculated. The doubling or halving continues until the 
step becomes too large to improve the result, or until it 
10 becomes small enough to not degrade the result. The resulting 
suitable step size is saved with d as part of the subgradient 
bundle. The last acceptable u + is assigned to u k +1 . 

Three distinct criteria are used to determine when u k is 
close enough to u k : 
15 1. The Wolfe's algorithm criterion of d=0 given that the 
test has been repeated with the bundle containing only 
the most recent subgradient . 

2. The difference between the lower bound @ k {u K ) and the 
upper bound v^z*, u k ) being less than a preset relative 

2 0 threshold. (Six percent was found to be an effective 

threshold for radar tracking problems.) 

3. An iteration count being exceeded. 
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The use of limits on the iteration are particularly- 
effective for iterations at the level 3 through (n-1) as these 
iterations will be repeated so as to resolve higher order 
coefficient sets. With limited iterations the process is in 
5 general robust enough to improve the estimate of upper order 
Lagrangian Coefficients. By limiting the iteration counts 
then the total processing time for the algorithm becomes 
deterministic. That characteristic means the process can be 
effective for real time problems such as radar, where the 
10 results of the last scan of the environment must be processed 
prior to the next scan being received. 

1-2.4. Solution Recovery 
The following process determines if the u k values at what 

15 is believed to be the peak of the function W k adequately 
approximate the u k values at the peak of the corresponding 
function & k . This is done by determining if the corresponding 
penalized cost function is minimized. Thus, the u k values at 
the peak of the function W k are first substituted into Problem 

20 Formulation (1.5) to determine a set of z assignments for k-1 
points. During the foregoing reduction process, Problem 
Formulation (1.7) yielded a tentative list of k selected z 
points that can be described by their indices as: 
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\ . . ijt-ij| W ° i / where N 0 is the number of cost elements selected 
into the solution. One possible solution of Problem 
Formulation (1.4) is the solution of Problem Formulation (1.5) 

^0 



with mf as it 



1 J = l 

was defined in Problem Formulation (1.6). If this solution 



satisfies the k th constraint set, then it is the optimal 
solution. 

However, if the solution for Problem Formulation (1.5) is 
not feasible (decision 355) , then the following adjustment 
10 process determines if a solution exists which satisfies the Jt^ 
constraint while retaining the assignments made in solving 
Problem Formulation (1.7) . To do this, a two-dimensional cost 
matrix is defined based upon all observations from the set 
which could be used to extend the relaxed solution. 

15 

k for 1 = 0, . . . f N k 

h j r c i 1/ ^. f i k _ 1/ i and j = 0f m m mfN ^ m (1.12) 

If the resulting two-dimensional assignment problem, 
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Minimize 
Subject to: 



5 

has a feasible solution, then the indices of that solution map 
to the solution of Problem Formulation (1.4) for the k- 
dimensional problem. The first index in each resultant 
solution entry is the pointer back to an element of the 

10 |^ il ? . . . ijt-i^^ ±k jj list. That element supplies the first Jc-1 
indices of the solution. The second index of the solution to 
the recovery problem is the k th index of the solution. 
Together these indexes specify the values of z k that solve 
Problem Formulation (1.4) at the k th level. 

15 If Problem Formulation (1.13) does not have a feasible 

solution, then the value of which was thought to represent 
Up is not representative of the actual peak and further 
iteration at the k th level is required. This decision 
represent the branch path from step 214 and equivalent steps. 

20 



j=Q 1=0 J 
1 = 0 J 



j=l, . . . ,N 0I 



w e {0,1} j= 0/ . 



,N 0 , 

.,N 0 , 1 = 0, 



(1.13) 
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1.3. Partitioning 
A partitioning process is used to divide the cost matrix 
that results from the Scoring step 154 into as many 
independent problems as possible prior to beginning the 
5 relaxation solution. The partitioning process is included 
with the Problem Formulation Step 310. The result of 
partitioning is a set of problems to be solved, i.e., there 
will be p 1 problems that consist of a single hypothesis, p 2 
problems that consist of two hypothesis, etc. Each such 

10 problem is a group in that one or more observations or tracks 
are shared between members of the group. 

In partitioning to groups no consideration is given to 
the actual cost values. The analysis depends strictly on the 
basis of two or more cost elements sharing the same specific 

15 axis of the cost matrix. In a two-dimensional case two cost 
elements must be in the same group if they share a row or a 
column. If the two elements are in the same row, then each 
other cost element that is also in the same row, as well as 
any cost elements that are in columns occupied by members of 

20 the row must be included in the group. The process continues 
recursively. In literature it is referred to as "Constructing 
a Spanning Forest." The 7c-dimensional case is analogous to 
the two-dimensional case. The specific method we have 
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incorporated is a depth first search, presented by Aho, 
Hopecraft, and Ullman (A.V. Aho, J.E. Hopcroft, and J.D. 
Ullman. Design and Analysis of Computer Algorithms . Addison - 
Wesley, MA, 1974) . 
5 The result of partitioning at level M is the set of 

problems described as {P i:f | i=l, . . . t Pj and j=l, . . . where N is 

the total number of hypothesis. The problems labeled 
{P il \i=l f ...,p 1 } are the cases where there is only one choice 
for the next observation at each scan and that observation 

10 could be used for no other track, i.e., it is a single 
isolated track. The hypothesis must be included in the 
solution set and no further processing is required. 

As hypothesis are constructed the first element is used 
to refer to the track ID. Any problem in the partitioned set 

15 which does not have shared costs in the first scan represent 
a case where a track could be extended in several ways but 
none of the ways share an observation with any other track. 
The required solution hypothesis for this case is the 
particular hypothesis with the maximum likelihood. For this 

20 case all likelihoods are determined as was described in 
Scoring and the maximum is selected. 

In addition to partitioning at the level M, partitioning 
is applied to each subsequent level M-l,...,2. For each 
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problem that was not eliminated by the prior selection, 
partitioning is repeated, ignoring observations that are 
shared in the last set of observations. Partitioning 
recognizes that relaxation will eliminate the last constraint 
5 set and thus partitioning is feasible for the lower level 
problems that will result from relaxation. This process is 
repeated for all levels down to k=3 . The full set of 
partitionings can be performed in the Problem Formulation Step 
310, prior to initiating the actual relaxation steps. The 

10 actual two-dimensional solver used in step 206 includes an 
equivalent process so no advantage would be gained by 
partitioning at the k=2 level . 

There are two possible solution methods for the remaining 
problems. "Branch and Bound" as was previously described, or 

15 the relaxation method that this invention describes. If any 
of the partitions have 5-10 possible tracks and less than 50 
to 2 0 hypotheses, then the prior art "Branch and Bound" 
algorithm generally executes faster than does the relaxation 
due to its reduced level of startup overhead. The "Branch and 

2 0 Bound" algorithm is executed against all remaining M level 
problems that satisfy the size constraint. For the remainder 
the Relaxation algorithm is used. The scheduling done in 
Problem Formulation allows each Problem Formulation (1.5) cost 
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matrix resulting from the first step of relaxation to be 
partitioned. The resulting partitions can be solved by any of 
the four methods: isolated track direct inclusion, isolated 
track tree evaluation, small group "Branch and Bound" or an 
5 additional stage of relaxation as has been fully described. 

The partitioning after each level of Lagrangian 
Relaxation is effective because when a problem is relaxed, the 
constraint that requires each point to be assigned to only one 
track is eliminated (for one image at a time) . Therefore, two 

10 tracks previously linked by contending for one point will be 
unlinked by the relaxation which permits the point to be 
assigned to both tracks. The fruitfulness of partitioning 
increases for lower and lower levels of relaxation. 

The following is a more detailed description of the 

15 partitioning method. Its application at all but the level M 
depends upon the relaxation process described in this 
invention. The recursive partitioning is therefore a distinct 
part of this invention. The advantage of this method is 
greatly enhanced by the sparse cost matrix resulting from 

2 0 tracking problems. However the sparse nature of the problem 
requires special storage and search techniques. 

A hypothesis is the combination of the cost element c n the 
selection variable z n and all observations that made up the 
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potential track extension. It can be written as, 

h a ={c n ,z n , {o n] =o kl \k=l f . . . ,M ± e {1, . . . ,N k }}} , i.e., cost, 
selection variable and observations in the hypothesis. Here 
n . . . f N} , where N is the total number of hypotheses in the 

5 problem. While the cost and assignment matrices were 
previously referenced, these matrices are not effective 
storage mechanisms for tracking applications. Instead the 
list of all hypothesis and sets of lists for each dimension 
that reference the hypothesis set are stored. The 
10 hypothetical set in list form is: 



For each dimension k=l,...,M there exists a set of lists, 
with each list element being a pointer to a particular 
15 hypothesis: 

i=N k ,k=M 



where N k . is a number of hypothesis containing the i th 
observation from scan k. This structure is a multiply linked 
list in that any observation is associated with a set of 
pointer to all hypothesis it participates in, and any 
2 0 hypothesis has a set of pointers to all observations that 
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formed it. (These pointers can be implemented as true 
pointers or indices depending upon the particular programming 
language utilized. ) 

Given this description of the matrix storage technique 
5 then the partitioning technique is as follows: Mark the first 
list L k and follow out all pointers in that list to the 
indicated hypothesis h pk . for i=l, . . . , U k . Mark all located 
hypothesis, and for each follow pointers back the particular 
L k for Jc=l, . . . , AT. Those L f s if not previously marked get 

10 marked and also followed out to hypothesis elements and back 
to L's. When all such L's or h ] s being located are marked, 
then an isolated sub-problem has been identified. All marked 
elements can be removed from the lists and stored as a unique 
problem to be solved. The partitioning problem then continues 

15 by again starting at the first residual set L. When none 
remain, the original problem is completely partitioned. 

Isolated problems can result from one source track having 
multiple possible extensions or from a set of source tracks 
contending for some observations. Because one of the indices 

2 0 of k (in our implementation it is k=l) indicates the source 
track then it is possible to categorize the two problem types 
by observing if the isolated problem includes more than one 
L-list from the index level associated with tracks. 
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1.4. Subsequent Track Prediction 
As noted above, after the points are assigned to the 
respective tracks, some action is taken such as controlling 
take-offs and landings to avoid collision, advising airborne 
5 aircraft to change course, warning airborne aircraft of an 
adjacent aircraft in a commercial environment, or aiming and 
firing an anti-aircraft (or anti-missile) missile, rocket or 
projectile, or taking evasive action in a military 
environment. Also, the tracking can be used to position a 

10 robot to work on an object. For some or all of these 
applications, usually the tracks which have just been 
identified are extended or extrapolated to predict a 
subsequent position of the aircraft or other object. The 
extrapolation can be done in a variety of prior art ways; for 

15 example, straight line extensions of the most likely track 
extension hypothesis, parametric quadratic extension of the x 
versus time, y versus time, etc. functions that are the result 
of the filtering process described earlier, least square path 
fitting or Kalman filtering of just the selected hypothesis. 

20 The process of gating, filtering, and gate generation as 

they were previously described require that a curve be fit 
through observations such that fit of observations to the 
likely path can be scored and further so that the hypothetical 
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target future location can be calculated for the next stage of 
gating. In the implementation existing, quadratic functions 
have been fit through the measurements that occur within the 
window. A set of quadratics is used, one per sensor 
5 measurement. To calculate intercept locations these 

quadratics can be converted directly into path functions. 
Intercept times are then calculated by prior methods based 
upon simultaneous solution of path equations. 

The use of the fitted quadratics is not as precise as 

10 more conventional filters like the Kalman Filter. They are 
however much faster and sufficiently accurate for the relative 
scores required for the Assignment problem. When better 
location prediction is required then the assignment problem is 
executed to select the solution hypothesis and based upon the 

15 observations in those hypothesis the more extensive Kalman 
filter is executed. The result is tremendous computation 
savings when compared with the Kalman Filter being run on all 
hypothetical tracks . 

Based on the foregoing, apparatus and methods have been 

2 0 disclosed for tracking objects. However, numerous 

modifications and substitutions can be made without deviating 
from the scope of the present invention . For example, the 
foregoing functions ¥ can also be defined as recursive 
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approximation problem in which several values of higher order 
u k values are used to eliminate the higher than approximation 
characteristic of the function <P k . The hill climbing of the 
function ¥ can be implemented by a high order hill climbing 
5 using the enhanced function & k . Although the result would not 
be as efficient it seems likely that the method would 
converge. Therefore, the invention has been disclosed by way 
of example and not limitation, and reference should be made to 
the following claims to determine the scope of the present 
10 invention. 



II. AN ALTERNATIVE EMBODIMENT OF THE MULT I -DIMENSIONAL 
ASSIGNMENT SOLVING PROCESS 

The following description provides an alternative second 

15 embodiment of the multi-dimensional assignment solving process 

of section 1 . 1 . However, in order to clearly describe the 

present alternative embodiment, further discussion is first 

presented regarding the data assignment problem of 

partitioning measurements into tracks and false alarms and, in 

20 particular, the representation of multi-dimensional assignment 

problems . 
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II. 1. Formulation of the Assignment Problem 
The goal of this section is to explain the formulation of 
the data association problems, and more particularly multi- 
dimensional assignment problems, that govern large classes of 
5 problems in centralized or hybrid centralized-sensor level 
multisensor/multitarget tracking. The presentation is brief; 
technical details are presented for both track initiation and 
maintenance in (A.B. Poore. Multidimensional assignment 
formulation of data association problems arising from 

10 multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications, 3:27-57, 1994) 
for non- maneuvering targets and in (A.B. Poore. 
Multidimensional assignments and multitarget tracking: 
Partitioning data sets. In P. Hansen, I.J. Cox, and B. 

15 Julesz, editors, DIMACS Series in Discrete Mathematics and 
Theoretical Computer Science, volume 19, pages 169-198, 
Providence, R.I., 1995. American Mathematical Society) for 
maneuvering targets. Note that the present formulation can 
also be modified to include target features (e.g., size and 

20 type) into the scoring process 154. 

The data assignment problems for multisensor and 
multitarget tracking considered here are generally posed as 
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that of maximizing the posterior probability of the 
surveillance region (given the data) according to 



Maximize i 



P[T= Y \Z-) . 



p{v=y°\z M ) 



(2.1) 



where Z M represents M data sets, y is a partition of indices 
of the data (and thus induces a partition of the data) , r* is 
the finite collection of all such partitions, r is a discrete 
random element defined on r* , y° is a reference partition, and 
10 P(F=y\z M ) is the posterior probability of a partition y being 
true given the data Z M . The term partition is defined below. 

Consider M observation sets Z{k) (k=l,... t N) each of N k 
observations \z t J\ , and let Z M denote the cumulative data set 
defined by 



Z(k) =|zi* J W * and Z m = {z(l) , . . . r Z(M)} f {2.2\ 



15 respectively. In multisensor data fusion and multitarget 
tracking the data sets Z{k) may represent different classes of 
objects, and each data set can arise from different sensors. 
For track initiation the objects are measurements that must be 
partitioned into tracks and false alarms. In the formulation 

2 0 of track extensions a moving window over time of observations 
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sets is used. The observation sets will be measurements which 
are: (a) assigned to existing tracks, (b) designated as false 
measurements, or (c) used for initiating tracks. However, 
note that alternative data objects instead of observation sets 
5 may also be fused such as in sensor level tracking wherein 
each sensor forms tracks from its own measurements and then 
the tracks from the sensors are fused in a central location. 
Note that, as one skilled in the art will appreciate, both 
embodiments of the present invention may be used for this type 

10 of data fusion. 

The data assignment problem considered presently is 
represented as a case of set partitioning defined in the 
following way. First, for notational convenience in 

representing tracks, a zero index is added to each of the 

15 index sets in (2.2), a gap filler z£ is added to each of the 
data sets Z(k) in (2.2), and a "track of data" is defined as 



partition of the data will refer to a collection of tracks of 
data or track extensions wherein each observation occurs 
2 0 exactly once in one of the tracks of data and such that all 
data is used up. Note that the occurrence of the gap filler 
is unrestricted. The gap filler serves several purposes in 
the representation of missing data, false observations, 




where i k can now assume the values of 0 . A 
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initiating tracks, and terminating tracks. Note that a 
"reference partition" may be defined which is a partition in 
which all observations are declared to be false. 

Next under appropriate independence assumptions it can be 
shown that : 

_^E^U Ls EL.., (2.3) 



where L^...^ is a likelihood ratio containing probabilities for 
detection, maneuvers, and termination as well as probability 
density functions for measurement errors, track initiation and 
termination. Then if c£ 2 ... lff = - In (L± lm . ,± M ) , 



-In 



P{r=y\Z M ) 



P(r = Y°|Z 



M 



(i lf ...,i M )€y 



C. 



(2.4) 



Using (2.3) and the zero-one variable , =1, if 

(i I7 ...,i M ) E y, and 0, otherwise, one can then write the 
problem (2.1) as the following M-dimensional assignment 
problem (as also presented in section Problem Formulation 
(1.4) with k=M) : 
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(a) Minimize 



i^O 2„=0 Ml M 



(b) Subject to: X - - . X ^ . =1/ i 2 = 1, . . . 



(c) 



£ • • • S £ • • • £ z i i =i > ( 2 - 5 ) 

for = 2, . . . ,N k and k=2, . . . ,M-1, 

x v • • • / ^ Z> ■ 1^ — 1/ . . . fN Mf 

™ . _ • • • -l-ir 



i =0 i =0 1 " M 

-4 u z /f-i u 



5 (d) 

(e) Z V-.^ 6 {0 ' l} for all i x ,...,i M , 

where c 0 ... 0 is arbitrarily defined to be zero. Here, each 
group of sums in the constraints (2.5) (b)-(2.5) (e) represents 

10 the fact that each non-gap filler observation occurs exactly 
once in a "track of data." One can modify this formulation to 
include multi -assignments of one, some, or all the actual 
observations. The assignment problem (2.5) is changed 
accordingly. For example, if z£ is to be assigned no more 

15 than, exactly, or no less than n^ k times, then the "=l" is 
changed to "<, = , > n k ±kt " respectively. 

Expressions for the likelihood ratios I/f 2 ... iM appearing in 
the costs c. . =-ln(L. . ) are well-known. In particular, 
discussions of these ratios may be found in (A . B . Poore . 

20 Multidimensional assignment formulation of data association 
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problems arising from multitarget tracking and multisensor 
data fusion, Computational Optimization and Applications, 
3: 27-57 , 1994; A.B. Poore . Multidimensional assignments and 
multitarget tracking: Partitioning data sets. In P. Hansen, 
5 I.J. Cox, B. Julesz, editors, DIMACS Series in Discrete 
Mathematics and Theoretical Computer Science, volume 19, pages 
169-198, Providence, R.I., 1995. American Mathematical 
Society) . Furthermore, the likelihood ratios are easily 
modified to include target features and to account for 

10 different sensor types. Also note that in track initiation, 
the M observation sets provide observations from M sensors, 
possibly all the same. Additionally note that for track 
maintenance, a sliding window of M observation sets and one 
data set containing established tracks may be used. In this 

15 latter case, the formulation is the same as above except that 
the dimension of the assignment problem is now M+l. 



II. 2. Overview of the New Lagrangian 
Relaxation Algorithms 

2 0 Having formulated an M-dimensional assignment problem 

(2.5), we now turn to a description of the Lagrangian 

relaxation algorithm. Subsection II. 2.1 below presents many 

of the relaxation properties associated with the relaxation of 

an n-dimensional assignment problem to an in-dimensional one 
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via a Lagrangian relaxation of n-m sets of constraints wherein 
m < n < M and preferably in the present invention embodiment 
n-m >1. Although any n-m constraint sets can be relaxed, the 
description here is based on relaxing the last n-m sets of 
5 constraints and keeping the first m sets. Given either an 
optimal or suboptimal solution of the relaxed problem, a 
technique for recovering a feasible solution of the 
n-dimensional problem is presented in subsection II. 2. 2 below 
and an overview of the Lagrangian relaxation algorithm is 
10 given in subsection II. 2. 3. 

The following notation will be used throughout the 
remainder of the work. Let M be an integer such that M >3 and 
let n£{3 , . . . , M} . The n-dimensional assignment problem is 



(a) Minimize: v n (z) 




15 



(b) Subject To 




t • • • 



(2.6) 



(c) 




for i k = l 



f • 



• r 



N k and k = 2 



f * * • r 



n-l r 



20 



(d) 




(e) 



z^.^ ejo,!} for all i 



1' • ' • ' i n- 
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To ensure that a feasible solution of (2.6) always exists, all 
variables with exactly one nonzero index (i.e., variables of 
the form zg... oio ... 0 for i k ±0) are assumed free to be assigned 
and the corresponding cost coefficients are well-defined. 



II. 2.1. The Lagrancrian Relaxed Assignment Problem 
The n-dimensional assignment problem (2.6) has n sets of 
constraints. A (N k +1) -dimensional multiplier vector associated 
with the k th constraint set will be denoted by u k -{u^ f u\ r . . . , u^ k ) T 

10 with Uq=0 and k=l,...,n. The n-dimensional assignment problem 
(2.6) is relaxed to an zn-dimensional assignment problem by 
incorporating n-m of the n sets of constraints into the 
objective function (2.6) (a) . Although any constraint set can 
be relaxed, the description of the relaxation procedure for 

15 (2.6) will be based on the relaxation of the last n-m sets of 
constraints. The relaxed problem is 



5 



u n ) = Minimize d> (z n ;u 



(2.7) 
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2 1 =0 i n =0 L k=m+l J - 

X • ■ • 22 z l...± n =1 ' i^l, 



Subject To 



for 1^ = 1, . . . , N k and k= 2, . . . ,m, 
z i 1 ...i n E { Q f 1 } for a11 i-i* • • • /i n / 



wherein we have added the multiplier Uq=0 for notational 
convenience. Thus, the multiplier u*e]R N * +1 with u£=0 is fixed. 

An optimal (or suboptimal) solution of (2.7) can be 
constructed from that of an in- dimensional assignment problem. 
To show this, define for each . . . , i m ) an index (j m+1/ • . * ,j n ) 

Um+i (^i/ * • • / in) / • ■ • / 3n • • * / in) ) &nd a new cost function 

* * * 'J J = 

-argmin <c* i + £ I i , = 0, . . . , N, and k = m + l, . . . ,nL (2.8) 
c i 1 ...i.= c i ^ , + 2^ Uj t for (1 ,...,! )* (0,...0), 
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c£..o = £ • • • £ Minimum jo,c 0 n _ 0i ± n + £ u£ k 

If (jjn+i/ • • • / Jn) is not unique, choose the smallest such j m+1 , 
amongst those {n-m) -tuples with the same j m+1 choose the 
smallest j m+2 , etc., so that ( j m+1/ . . . , j n ) is uniquely defined.) 
Using the cost coefficients defined in this way, the following 
in-dimensional assignment problem is obtained: 
0 mn (u m+1 / . . . f u n ) = Minimize 0^ (z^/u^ 1 , . . . , u n ) = v m (z m ) 

in m 



L=0 i m =0 
l m 



Subject to: 



N 0 N 

E • • • £ <...i m =1 ' ii^l/.-./Wi, (2.9) 



i 2 =0 i„=0 



i 1 = 0 i H -0 i k+1 = 0 i a =Q 

10 for = 1, . . . ,2^ and Jc=2, . . . / m-l / 

• * " l—J ^ . . . i — I ' 3-m~-i /•••/ -^m/ 

i,=0 i ,=0 1 

l m-l 



<...i M e {0,1} for all i lf . . .,i m . 



As an aside, observe that any feasible solution z n of (2.6) 
15 yields a feasible solution z m of (2.9) via the construction 



z m ■ = < 



1/ if ^...Vxin = 1 f ° r SOme ( W 

0, otherwise. 
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Thus, the m- dimensional assignment problem (2.9) has at least 
as many feasible solutions of the constraints as the original 
problem (2.6) . 

The following Fact A. 2 has been shown to be true. It 
5 states that an optimal solution of Problem Formulation (2.7) 
can be computed from that of Problem Formulation (2.9). 
Furthermore, the converse to this fact is provided in Fact 
A. 3. Moreover, if the solution of either of these two 
problems (2.7) or (2.9) is E-optimal (i.e., the objective 
10 function associated with the suboptimal solution is within "G" 
of the objective function associated with the optimal 
solution) , then so is the other. 

Fact A. 2. Let W" be a feasible solution to problem (2.9) 

and define vf 1 by 

15 w^...^*^...^ if (i m+1 , . . -i n ) = (j m+ i, . . .jj and 
. . . , ij * (0, . . . , 0) 

^...^=0 if (i n+1 , . . .i a ) * (j m+1 , . . .j n ) and (2.10) 

(i 2/ . . . ,ij # (0, . . . , 0) 

20 wg... ow ... lfl =l if c3... 0W ... Jb + £ < <0 

oi , i=0 if cj oi , , + /J >0. 

Then w 12 is a feasible solution of the Lagrangian relaxed 
problem (2.7) and 
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,u°) = <P m AWm** 1 , 



,u n >- £ £ u* 



k=m+l i k =0 



If, in addition, is optimal for (2.9), then is an optimal 
solution of (2.7) and 

^(u-v..,u*) = <p an ( U - +i / ...,u a )- £ i; < . 



5 With the exception of one equality being converted to an 

inequality, the following Fact is a converse of Fact A. 2 and 
has also been shown to be true. 

Fact A. 3. Let / be a feasible solution to problem (2.7) 
and define vf 1 by 

10 Wi^.V* ll w? 2 _ in for (i 2/ . . .,i m ) * (0, . . . ,0) , 

wS... 0 =0 if (i 2/ . . .,ij = (0, . . .,0) and 



...oi ffl+1 ...i n + £ "i* >° for all (i m+1 , . . , (2.11) 



Jt=JH+l ^ 



wS... 0 =l if (i 2/ • . . ,i m ) = (0, . . . , 0) and 
^-.-oi^.-.in + £ u i ^° for some (i^i/ • • -'i n > - 



15 
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Then vf 1 is a feasible solution of the problem (2.9) and 
0 ma (tf;vr 1 ,...,u a ) > <f mn (w m ;u m+1 r ... f u n )-^ £ < • 



If, in addition, vf 1 is optimal for (2.7) , then w" 1 is an optimal 
solution of (2.9) , 

<p mn (w n ;u m+1 ,...,u n ) =<p mn (w m ;u* +l ,...,u n ) - £ £ u k i r and 

k=m+l i k =0 

0 [u m+1 f . . .,u n ) =<P (u m+1 , . . .,u n ) " £ £ Ui*. 

5 II. 2. 2. The Recovery Procedure 

The next objective is to explain a recovery procedure or 
method for recovering a solution to the n-dimensional problem 
of (2.6) from a relaxed problem having potentially 
substantially fewer dimensions than (2.6). Note that this 

10 aspect of the alternative embodiment of the present invention 
is substantially different from the method disclosed in the 
first embodiment of the multidimensional assignment solving 
process of section 1.1 of this specification. Further, this 
alternative embodiment provides substantial benefits in terms 

15 of computational efficiency and accuracy over the first 
embodiment/ as will be discussed. Thus, given a feasible 
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(optimal or suboptimal) solution vf" of (2.9) (or vf 1 if Problem 
Formulation (2.7) is constructed via Fact A. 2) , an explanation 
is provided here regarding how to generate a feasible solution 
z n of (2.6) which is close to / in a sense to be specified. 
We first assume that no variables in (2.6) are preassigned to 
zero (this assumption will be removed shortly) . The 
difficulty with the solution is that it need not satisfy the 
last n-m sets of constraints in (2.6) . Note, however, that if 
vf 1 is an optimal solution for (2.9) and w 21 (constructed as in 
Fact A. 2) satisfies the relaxed constraints, then vf 1 is optimal 
for (2.6) . The recovery procedure described here is designed 
to preserve the 0-1 character of the solution vf" of (2.9) as 
far as possible. That is, if w^... ij77 =l and i 2 *0 for at least one 
1=1, . . . f m, then the corresponding feasible solution z n of (2.6) 
is constructed so that tf± lmmmiD =l for some (i ffl+1/ . . . / i n ) . Note 
that by this reasoning, variables of the form z2... oijn+1 .„ lj2 can be 
assigned to a value of 1 in the recovery problem only if 
w£.. 0 =l. However, variables &o...oi m+1 ...i n wil1 be treated 
differently in the recovery procedure in that they can be 
assigned 0 or 1 independent of the value wS... 0 . This increases 
the feasible set of the recovery problem, leading to a 
potentially better solution. 
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Let | (i{ t - . • , ±1) Y\ be an enumeration of indices of w m (or 
the first m indices of W 2 constructed as in Fact A. 2) such that 
w m j , = 1 and (il, . . . , i j a ) * (0, . . . ,0) . Set {i° lf . . . , i°> = (0 , . . . , 0) 
for j-0 and define 

5 cjT +1 , = c fl , ... (2.12) 

for 1^=0,...,^; Jc=m+1 , . . . ,n; j=0,...,N g . 

Let Y denote the solution of the (n-m+1) -dimensional 
10 assignment problem: 

Minimize 22 • ■ • z2 c ji m . . .± Yu i 
Subject to 2J ■ • • X) ,• = 1 ' J =1 ' • • • '^o/ (2.13) 

L L --L ••■2. y Jiw ...i n =1 ' 

for ijc=2, . . . ,N k and k=m+l, . . . ,n-l, 

E E • ' • L y« i =1 ' i n= 1 ' • • • ' N *> 

•^ji ntl ...i n e for all j, i mlf ...,i a . 



The recovered feasible solution z 11 of (2.6) corresponding to 
20 the multiplier set {u m , . . . ,u n } is then defined by 
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n 

z i,...i r 



1, if {i 1 ,...,i m ) = (ii,... f il) for some 

j = 0, . . .,N 0 and Y . ..i a =l, (2.14) 
0, otherwise. 



This recovery procedure is valid as long as all cost 
coefficients c n are defined and all zero-one variables in z n 
are free to be assigned. Note that modifications are 
necessary for sparse problems. If the cost coefficient 
c n j tJt .is well defined and the zero-one variable 



c n j tjt .is free to be assigned to zero or one, then define 
c n-m+i n as in {2 12) wit h z.: m+1 ± being free to 

be assigned. Otherwise, z 7 n r m+1 i is preassigned to zero. To 
ensure that a feasible solution exists, we now only need 
10 ensure that the variables ZjoTT.o are free for j=0, 1, . . . ,N 0 . 
(Recall that all variables of the form zS.. iijt ... 0 are free (to be 
assigned) and all coefficients of the form c?... ijc ... 0 are well 
defined for Jc=l, . . . ,n. ) This is accomplished as follows. If 
the cost coefficient 7 - 1 is well defined and 

15 ^^..i^.^o is fre6/ then d6f lne c J°*" 0 = C i 1 h 2 \..±lo...o With z "o!".?o 
being free. Otherwise, since all variables of the form 

Zo...i k ...o are known feasible and have well-defined costs, put 

n-m+l _ V A n 
C jC...0~ 0...0i, J 0...0 * 
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II. 3. The Mult i -Dimensional Lagrangian Relaxation 

Algorithm for the n-Dimensional Assignment 
Problem 

Starting with the M- dimensional assignment problem (2.6) , 
5 i.e., n=M, the algorithm described below is recursive in that 
the M-dimensional assignment problem is relaxed to an tri- 
dimensional problem by incorporating {n-m) sets of constrains 
into the objective function using Lagrangian relaxation of 
this set. This problem is maximized with respect to the 

10 Lagrange multipliers, and a good suboptimal solution to the 
original problem is recovered using an {n-m+1) -dimensional 
assignment problem. Each of these two (the m-dimensional and 
the (n-m+1) -dimensional assignment problems) can be solved in 
a similar manner. 

15 More precisely, reference is made to Figure 8 which 

presents a flowchart of a procedure embodying the multi- 
dimensional relaxation algorithm, referred to immediately 
above. This procedure, denoted MULT I _D I M_RE LAX in Figure 8, 
has three formal parameters , n, PROB_FORMULAT I ON and 

20 ASSIGNMT_SOLUTION, which are used to transfer information 
between recursive instantiations of the procedure. In 
particular, the parameter, n, is an input parameter supplying 
the dimension for the mult i -dimensional assignment problem (as 
in (2. 6)) which is to be solved (i.e., to obtain an optimal or 
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near-optimal solution) . The parameter, PROB_FORMULAT I ON , is 
also an input parameter supplying the data structure (s) used 
to represent the n-dimensional assignment problem to be 
solved. Note that PROB___FORMULAT I ON at least provides access 
5 to the cost matrix, [c n ] , and the observation sets whose 
observations are to be assigned. Additionally, the parameter, 
ASSIGNMT_SOLUTION, is used as an output parameter for 
returning a solution of a lower dimensioned assignment problem 
to an instantiation of MULTI_DIM_RELAX which is processing a 

10 higher dimensioned assignment problem. 

A description of Figure 8 follows. Assuming 
MULT I _D I M_RE LAX is initially invoked with M as the value for 
the parameter n and the PROB_FORMULAT I ON having a data 
structure (s) representing an M-dimensional assignment problem 

15 as in (2.5), in step 500 an integer m, 2<m<n is chosen such 
that the constraint sets corresponding to dimensions . . . ,n 

are to be relaxed so that an m-dimensional problem formulation 
as in (2.7) results. In step 504 an initial approximation is 
provided for { u£ +1 , • . . , u n 0 ) . Subsequently, in step 508 the above 

20 initial values for { u^ 1 , . . . , u"} are used in an iterative 
process in determining (£P +I , ...,£?) which maximizes 

{^(u™* 1 , . . . ,u n ) where u k GW^ +1 and k=m+l, . . . ,n} for a feasible 
solution w 13 subject to the constraints of Problem Formulation 
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(2.7). Note that by maximizing ^ {u m+1 , . . . , u n ) , some of the 
constraints that were relaxed are being forced to be satisfied 
and in so doing information is built into a solution of this 
function for solving the input assignment problem in 
5 11 PROB_FORMULATION . " Further note that a non- smooth 

optimization technique is used here and that a preferred 
method of determining a maximum in step 508 is the bundle- 
trust method of Schramm and Zowe (H. Schramm and J. Zowe. A 
version of the bundle idea for minimizing a non-smooth 
10 function: Conceptual idea, convergence analysis, numerical 
results. SIAM Journal on Optimization, 2, No. 1:121-152, 
February, 1992) . This method, along with various other 
methods for determining the maximum in step 508, are discussed 
below. 

15 Also note that for m >2 , a solution to the optimization 

problem of step 508 is NP-hard and therefore cannot be solved 
optimally. That is, there is no known computationally 
tractable method for guaranteeing an optimal solution. Thus, 
there are two possibilities: either (a) allow m to be greater 

2 0 than 2 and use auxiliary functions similar to those disclosed 
in the first embodiment of the Jc-dimensional assignment solver 
300 in section I to compute a near-optimal solution, or (b) 
make m=2 wherein algorithms such as the forward/ reverse 
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auction algorithm of D.P. Bertsekas and D.A. Castanon (D.P. 
Bertsekas and D.A. Castanon. A forward/ reverse auction 
algorithm for asymmetric assignment problems. Computational 
Optimization and Applications, 1: 277-298, 1992) provides an 
5 optimal solution. 

If option (a) immediately above is chosen, then the 
auxiliary functions are used as approximation functions for 
obtaining at least a near-optimal solution to the optimization 
problem of step 508. Note that the auxiliary functions used 

10 depend on the value of m. Thus, auxiliary functions used when 
m=3 will likely be different from those for m=4 . But in any 
case, the optimization procedure is guided by using the merit 
function, 0 2n (u 3 , . . . , u n ) , which can be computed exactly via a 
two-dimensional assignment problem for guiding the 

15 maximization process. 

Alternatively, if option (b) above is chosen, then two 
important advantages result, namely, the optimization problem 
of step 508 can be always solved optimally and without using 
auxiliary approximation functions. Thus, better solutions to 

2 0 the original M-dimensional assignment problem are likely since 
there is no guarantee that when the non- smooth optimization 
techniques are applied to the auxiliary functions the 
techniques will yield an optimal solution to step 508. 



95 



CSUROWSRf 

Furthermore, it is important to note that without auxiliary- 
functions, the processing in step 508 is both conceptually 
easier to follow and more efficient. 

Subsequently, in step 512 of Fig. 8, the solution 
5 (i?* 1 , . . . , u n ) is used in determining an optimal solution w" 1 to 
Problem Formulation (2.9) as generated according to (2.8) and 
Fact A. 3 . 

In step 516, the solution w 12 is used in defining the cost 
matrix c n ~ m+1 as in (2.12). Subsequently, if n-m+l=2, then the 

10 assignment problem (2.13) may be solved straightforwardly 
using known techniques such as forward/reverse auction 
algorithms. Following this, in step 528, the solution to the 
two-dimensional assignment problem is assigned to the variable 
ASSIGNMT_SOLUTION and in step 532 ASSIGNMT_SOLUTION is 

15 returned to a dimension three level recursion of 
MULTI_DIM_RELAX for solving a three-dimensional assignment 
problem. 

Alternatively, if in step 520, n-m+l>2 , then in step 536 
the data structure (s) representing a problem formulation as in 
20 (2.13) is generated and assigned to the parameter, 
PROB_FORMULATION. Subsequently, in step 54 0 a recursive copy 
of MUL T I _D I M_R E LAX is invoked to solve the lower dimensional 
assignment problem represented by PROB_FORMULATION. Upon the 
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completion of step 540, the parameter, ASSIGNMT_SOLUTION, 
contains the solution to the {n-m+1) -dimensional assignment 
problem. Thus, in step 544, the {n-m+1) th solution is used to 
solve the n-dimensional assignment problem as discussed 
5 regarding equations (2.14) . Finally, in steps 548 and 552 the 
solution to the n-dimensional assignment problem is returned 
to the calling program so that, for example, it may be used in 
taking one or more actions such as (a) sending a warning to 
aircraft or sea facility; (b) controlling air traffic; (c) 
10 controlling anti-aircraft or anti-missile equipment; (d) 
taking evasive action; or (e) surveilling an object. 



II. 3.1. Comments on the Various Procedures 
Provided by Figure 8 

15 

There are many procedures described by Figure 8 . One 
such procedure is the first embodiment of the multidimensional 
assignment solving process of section 1.1. That is, by 
specifying m-n-1 in step 500, a single set of constraints is 

20 relaxed in step 508. Thus, one set of constraints is 
incorporated into the objective function via the Lagrangian 
problem formulation, resulting in an (m=n-l) -dimensional 
problem. The relaxed problem is subsequently maximized in 
step 512 with respect to the corresponding Lagrange 

25 multipliers and then a feasible solution is reconstructed for 
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the n-dimensional problem using a two-dimensional assignment 
problem. The second procedure provided by Figure 8 is a novel 
approach which is not suggested by the first embodiment of 
section 1.1. In fact, the second procedure is somewhat of a 
5 mirror image of the first embodiment in that n-2 sets of 
constraints are simultaneously relaxed, yielding immediately 
an (m=2) -dimensional problem in step 512. Thus, a feasible 
solution to the n -dimensional problem is then recovered using 
a recursively obtained solution to an (n-1) -dimensional 

10 problem via step 540. In this case, the function values and 
subgradients of @ 2n (u 3 , . . . , u n ) of step 508 can be computed 
optimally via a two-dimensional assignment problem. The 
significant advantage here is that there is no need for the 
merit or auxiliary functions, W kf as required in the first 

15 embodiment of the multidimensional assignment solving process 
of section 1.1 above and also there is no need for the more 
general merit or auxiliary functions W^, as discussed in 
subsection II. 4. 2 below. Further, note that all function 
values and subgradients used in the nonsmooth maximization 

20 process are computed exactly (i.e., optimally) in this second 
procedure. Moreover, problem decomposition is now carried out 
for the n-dimensional problem; however, decomposition of the 
(n-1) -dimensional recovery problem (and all lower recovery 
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problems) is performed only after the problem is formulated. 

Between these two procedures are a host of different 
relaxation schemes based on relaxing n - m sets of constraints 
to an 27?-dimensional problem {2<m<n) , but these all have the 
5 same difficulties as the procedure for the first embodiment of 
section I . 1 in that the relaxed problem is an NP-hard problem. 
To resolve this difficulty, we use an auxiliary or merit 
function as described in subsection II. 4. 2 below. (The 

notation W k was used for W n _ lfI1 in section 1.1.) For the case 

10 m<n-l, the recovery procedure is still based on an NP-hard 

(jn-m+l) -dimensional assignment problem. Note that the 
partitioning techniques similar to those discussed in Section 
1.3 may be used to identify the assignment problem with a 
layered graph and then to find the disjoint components of this 

15 graph. In general , all relaxed problems can be decomposed 
prior to any nonsmooth computations because their structure 
stays fixed throughout the algorithm of Fig. 8. However, all 
recovery problems cannot be decomposed until they are 
formulated, as their structure changes as the solutions to the 

2 0 relaxed problems change. 
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II. 4. Details and Refinements Relating to the 

Flowchart of Figure 8 

Further explanation is provided here on how various steps 

of Figure 8 are solved. Note that the refinements presented 

5 here can significantly increase the speed of the relaxation 

procedure of step 508. 

II. 4.1 Maximization of the Non-smooth Function 
® ^(u m+1 u°> 

10 One of the key steps in the Lagrangian relaxation 

algorithm in section II. 3 is the solution of the problem 

Maximize { ^ (u m+ \ . . . , u n ) : u k eW k+1 ; k= m+1 , . . . , n} , (2.15) 

15 where u£=0 for all Jc=m+1, . . . , n. The evaluation of 

@mn (u™ 1 , . . . , u n ) requires the optimal solution of the 
corresponding minimization problem (2.7). The following 
discussion provides some properties of these functions. 

Fact A. 4. Let u m+1 , . . . , u n be multiplier vectors associated 

20 with the {m+1) st through the n th set of constraints (2.6), let 
0^ be as defined in (2.7) , let V n (F 1 ) be the objective function 
value of the n-dimensional assignment problem in equation 
(2.6) , let z n be any feasible solution of (2.6) , and let z* 2 be 
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an optimal solution of (2.6). Then, ^ {u m+1 , . . . , u 11 ) is 
piecewise affine, concave and continuous in { u m+1 , . . . , u n } and 

5 (2.16) 

Furthermore , 

tfm-lndx"/^ 1 / • • ./U°) < ^(ll** 1 , . ..,12°) (2.17) 

for all u*ER M * +I with u*=0 and k=m, . . . ,n. 

10 

Most of the procedures for nonsmooth optimization are 
based on generalized gradients called subgradients , given by 
the following definition. 

Definition . At u= (u^ 1 , . . . , u n ) the set d0 mn {u) is called a 
15 subdifferential of 0^ and defined by 

5<P JBn (u)={g G R^ +1 x. . .xR^ +I |* m (w)-# m (u)^g r (w-u) (2.18) 

\/w E R^ +2 x. . .xR^ +1 } . 

20 A vector g E ^^(u) is called a subgradient. 

If z 23 is an optimal solution of (2.7) computed during 
evaluation of ^(u), differentiating with respect to u5 n 
yields the following i n th component of a subgradient g of ^(u) 
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< = £••• £ £ ...£*i...i.-i < 2 - 19 > 

for ijt=l, . . . f N k and k-m+1, . . . ,n. 

5 If z 23 is the unique optimal solution of (2.1), 3^ (u) ={gr} and 
^ is dif ferentiable at u. If the optimal solution is not 
unique, then there are finitely many such solutions, say 
^(1) , . . . ,^{K) . Given the corresponding subgradients , 
g 1 f ... f g K l the subdif f erential d<X>{u) is the convex hull of 
10 {g 1 f . . . ,3*} • 

II. 4. 2 Mathematical Description of 

a Merit or Auxiliary Function 

For real-time needs, one must address the fact that the 

15 non-smooth optimization problem of step 508 (Fig. 8) requires 

the solution of an NP-hard problem for m>2 . One approach to 

this problem is to use the following merit or auxiliary 

function to decide whether a function value has increased or 

decreased sufficiently in the line search or trust region 

2 0 methods: 
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m+1 




(2.20) 




if m = 2 

or (2.7) is solved optimally, 
. , u n ) , otherwise. 



where the multipliers ,i? that appear in lower order 

5 relaxations used to construct (suboptimal) solutions of the 
m-dimensional relaxed problem (2.7) have been explicitly 
included. Note that is well-defined since (2.9) can always 
be solved optimally if m=2 . For sufficiently small problems 
(2.7) or (2.9), one can more efficiently solve the NP-hard 

10 problem by branch and bound. This is the reason for the 
inclusion of the first case; otherwise, the relaxed function 
0 2n is used to guide the nonsmooth optimization phase. That 
the merit function provides a lower bound for the optimal 
solution follows directly from Fact A. 4 and the following 

15 fact. 

Fact A. 5. Given the definition of in (2.20) above, 
^({?, . . . , tF/u"" 1 , . . . ,u n ) <0 mn {u m+1 / . . . ,u n ) for all multipliers 



if 1 ,!! 1 



.m+1 



(2.21) 



20 
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The actual function value used in the optimization phase is 
W^; however, the subgradients are computed as in (2.19), but 
with the solution z i 2 ...i n being a suboptimal solution 
constructed from a relaxation procedure applied to the m- 
5 dimensional problem. Again, the use of these lower order 
relaxed problems is the reason for the inclusion of the 
multipliers t?, . . . , u m . 

To explain how the merit function is used, suppose there 
is a current multiplier set ( u^, . . . , u n old ) and it is desirable 

10 to update to a new multiplier set . . . , u* ew ) via 

(uS£/ • • • / v&J = (u£S, . . . , u n old ) + (Zlu- 2 , . . . ,412°) . Then 
^ ( u*oiar • • . / "Sid/ u old/ • • • / "Sid) is computed where ( iP old , . . . , ug ld ) is 
obtained during the relaxation process used to compute a high 
quality solution to the relaxed zn-dimensional assignment 

15 problem (2.7) at (u m+1 , . . . , u n ) = (u^, . . . , u n old ) . The decision to 
accept (u££i/ - - - , u£ e J is then based on ( u* Id/ . . . , u™ ld ; 

u m 0 * ld , u n old ) < ( iP new/ u* ew ; u£i, . • . , o r some other 

stopping criteria commonly used in line searches. Again, 
( i?ew/ • • • / u£ e J represents the multiplier set used in the lower 

2 0 level relaxation procedure to construct a high quality 
feasible solution to the m-dimensional relaxed problem (2.7) 
at {if* 1 , . . . ,11") = (u£e£/ • * • f^new) • The point is that each time 
one changes ( u m+1 , . . . , u n ) and uses the merit function 
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¥^[1?, . . . , if 1 ; u m+1 , . . . f u n ) for comparison purposes, one must 
generally change the lower level multipliers {x?,... f u m ). 

An illustration of this merit function for m=n-l is given 
in (A.B. Poore and N. Rijavec. Partitioning multiple data 
5 sets: multidimensional assignments and Lagrangian relaxation. 
In P.M. Pardalos and H. Wolkowicz, editors, Qqudratic 
assignment and related problems: DIMACS Series in Discrete 
Mathematics and Theoretical Computer Science, volume 16, 
pages 25-37, 1994) . 

10 

II. 4. 3 Non-smooth Optimization Methods 
By Fact A. 4 the function ^(u) is a continuous, piecewise 
affine, and concave, so that the negative of ^(u) is convex. 
Thus the problem of maximizing ^(u) is one of nonsmooth 
15 optimization. Since there is a large literature on such 
problems, only a brief discussion of the primary classes of 
methods for solving such problems is provided. A first class 
of methods, known as subgradient methods, are reviewed and 
analyzed in (N.Z. Shor. Minimization Methods for Non- 
20 Differentiable Functions, Springer-Verlag, New York, 1985) . 
Despite their relative simplicity, subgradient methods have 
some drawbacks that make them inappropriate for the tracking 
problem. They do not guarantee a descent direction at each 
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iteration, they lack a clear stopping criterion, and exhibit 
poor convergence (less than linear) . 

A more recent and sophisticated class of methods are the 
bundle methods; excellent developments are presented in (J.-B. 
5 Hiriart-Urruty and C. Lemarechal . Convex Analysis and 
Minimization Algorithms I and II. Springer-Verlag, Berlin, 
1993; K.C. Kiwiel . Methods of Descent for Non-Dif f erentiable 
Optimization. In A. Dold and B. Eckmann, editors, Lecture 
Notes in Mathematics, volume 1133, Berlin, 1985. Springer- 

10 Verlag) . Bundle methods retain a set (or bundle) of 
previously computed subgradients to determine the best 
possible descent direction at the current iteration. Because 
of the nonsmoothness of 3^, the resulting direction may not 
provide a "sufficient 11 decrease in 3^. In this case, bundle 

15 algorithms take a "null" step, wherein the bundle is enriched 
by a subgradient close to u*. As a result, bundle methods are 
non-ascent methods which satisfy the relaxed descent condition 
^(Ufc+i) <<Pinn(Ujk) if u k+i*u k . These methods have been shown to 
have good convergence properties. In particular, bundle 

2 0 method variants have been proven to converge in a finite 
number of steps for piecewise affine convex functionals in 
(K.C. Kiwiel. An aggregate subgradient method for non-smooth 
convex minimization. Mathematical Programming, 27:320-341, 
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1983; H. Schramm and J. Zowe. A version of the bundle idea for 
minimizing a non- smooth function: Conceptual idea, convergence 
analysis, numerical results. SI AM Journal on Optimization, 2, 
No. 1:121-152, February, 1992). 
5 All of the above non- smooth optimization methods require 

the computation of ^(u) and a g 6 d^iii) . These in turn 
require the computation of the relaxed cost coefficients 
(2.8). In both the first and second procedures discussed in 
section 1 1. 3.1, maximization of 0 2 n( u ) must be repeatedly 

10 evaluated. In the most efficient implementations presently 
known of these two procedures, it was found that at least 95% 
of the computational effort of the entire procedure is spent 
in the evaluation of the relaxed cost coefficients (2.8) as 
part of computing 0 2 n( u ) • Thus, generally a method should be 

15 chosen that makes as efficient use of the subgradients and 
function values as possible, even at the cost sophisticated 
manipulation of the subgradients. In evaluating three 
different bundle procedures: (a) the conjugate subgradient 
method of Wolfe used in section 1 . 1 of the first embodiment of 

2 0 the present invention; (b) the aggregate subgradient method of 
Kiwiel (K.C. Kiwiel . An aggregate subgradient method for non- 
smooth convex minimization. Mathematical Programming, 27:320- 
341, 1983) ; and (c) the bundle trust method of Schramm and 
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Zowe (H. Schramm and J. Zowe. A version of the bundle idea for 
minimizing a non- smooth function: Conceptual idea, convergence 
analysis, numerical results. SIAM Journal on Optimization , 2, 
No. 1:121-152, February, 1992), it was determined that for a 
5 fixed number of non-smooth iterations, say, ten, the bundle- 
trust method provides good quality solutions with the fewest 
number of function and subgradient evaluations of all the 
methods, and is therefore the preferred method. 

10 II. 4. 4 The Two-Dimensional Assignment Problem 

The forward/ reverse auction algorithm of Bertsekas and 

Castanon (D.P. Bertsekas and D.A. Castanon. A forward/reverse 

auction algorithm for asymmetric assignment problems. 

Computational Optimization and Applications , 1:277-298, 1992) 
15 is used to solve the many relaxed two-dimensional problems 

that occur in the course of execution. 

II. 4. 5. Initial Multipliers and Hot Starts 
The effective use of "hot starts" is fundamental for 
20 real-time applications. A good initial set of multipliers can 
significantly reduce the number of non- smooth iterations (and 
hence the number of ^ evaluations) required for a high 
quality recovered solution. Further, there are several ways 
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that multipliers can be reused. First, if an n -dimensional 
problem is relaxed to an m-dimensional problem, relaxation 
provides the multiplier set {u^ 1 , u^ 2 , . . . , u n } . These can be 
used as the initial multipliers for the (n-m+l) -dimensional 
5 recovery problem for n-m+l>2 . This approach has also worked 
well to reduce the number of non- smooth iterations during 
recovery. 

Further, for track maintenance, initial feasible 
solutions are generated as follows. When a new scan of 

10 information (a new observation set) arrives from a sensor, one 
can obtain an initial primal feasible solution by matching new 
reports to existing tracks via a two-dimensional assignment 
problem. This is known as the track-while-scan (TWS) 
approach. Thus, an initial primal solution exists and then we 

15 use the above relaxation procedure to make improvements to 
this TWS solution. Also for track maintenance, multipliers 
are available from a previously solved and closely related 
multidimensional assignment problems for all but the new 
observation set. 

20 

II. 4. 6. Local Search Methods 
Given a feasible solution of the multidimensional 
assignment problem, one can consider local search procedures 
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to improve this result, as described in (C.H. Papadimitriou 
and K. Steiglitz. Combinatorial Optimization: Algorithms and 
Complexity. Prentice-Hall , Inc., Englewood Cliffs, NJ, 1982; 
J. Pearl. Heuristics: Intelligent Search Strategies for 
5 Computer Problem Solving. Addison-Wesley , Reading, MA, 1984) . 
One method is the idea of k interchanges. Since for sparse 
problems only those active arcs that participate in the 
solution are stored, it is difficult to efficiently identify 
eligible variable exchange sets with some data structures for 

10 solving the assignment problem. However, a local search that 
may be more promising is to use the two-dimensional assignment 
solver in the following way. Given a feasible solution to the 
multidimensional assignment problem, the indices that 
correspond to active arcs in the solution are enumerated. 

15 Subsequently, one coordinate is freed to formulate a two- 
dimensional assignment problem with one index corresponding to 
the enumeration and the other to the freed coordinate, and 
then solve a two-dimensional assignment problem to improve the 
freed index position. 



20 



II. 4. 7. Problem Decomposition 
The procedures described thus far are all based on 
relaxation. Due to the sparsity of assignment problems, 
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however, frequently a decomposition of the problem into a 
collection of disjoint components can be done wherein each of 
the components can be solved independently. Due to the setup 
costs of Lagrangian relaxation, a branch and bound procedure 
5 is generally more efficient for small components, say ten to 
twenty feasible 0-1 variables (i.e., 2^...^)- Otherwise, the 
relaxation procedures described above are used. Perhaps the 
easiest way to view the decomposition method is to view the 
reports or measurements as a layered graph. A vertex is 

10 associated with each observation point, and an edge is allowed 
to connect two vertices only if the two observations belong to 
at least one feasible track of observations. Given this 
graph, the decomposition problem can then be posed as that of 
identifying the connected subcomponents of a graph which can 

15 be accomplished by constructing a spanning forest via a depth 
first search algorithm, as described in (A.V. Aho, J.E. 
Hopcroft, and J.D. Ullman. The Design and Analysis of Computer 
Algorithms. Addison-Wesley, MA, 1974). Additionally, 
decomposition of a different type might be based on the 

20 identification of bi- and tri-connected components of a graph 
and enumerating on the connections. Here is a technical 
explanation. Let 

2 = {zi 1 i 2 „.ijzi 1 i 2 ...i a is not preassigned to zero] 
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denote the set of assignable variables. Define an undirected 
graph G(N,A) where the set of nodes is 

N={z n in \n=l, . . . ,n; i n =l, . . . ,N m } 

and arcs, 

5 A={ (z5 n ,zj 2 ) \n*l, j n *0 r j 2 *0 and there exists z ili2 ... in eZ 

such that j n =i n and Ji=ii} - 
Note that the nodes corresponding to zero index have not 
been included in the above defined graph, since two variables 
that have only the zero index in common can be assigned 
10 independently. Connected components of the graph are then 
easily found by constructing a spanning forest via a depth 
first search. (A detailed algorithm can be found in the book 
by Aho 7 Hopcroft and Ullman cited above) . Furthermore, this 
procedure is used at each level in the relaxation, i.e., is 
15 applied to each assignment problem (1.4) for n=3,...,M. 

The original relaxation problem is decomposed first. All 
relaxed assignment problems can be decomposed a priori and all 
recovery problems can be decomposed only after they are 
formulated. Hence, in the n-to- (n-1) case, we have n-2 
2 0 relaxed problems that can all be decomposed initially, and the 
recovery problems that are not decomposed (since they are all 
two-dimensional) . In the n-to-2 case, we have only one 
relaxed problem that can be decomposed initially. This case 



112 



CSUR.Q1USRI 

yields n-3 recovery problems, which can be decomposed only 
after they are formulated. 



Ill . NEW APPROACHES TO TRACK INITIATION AND MAINTENANCE 
5 USING MULTIDIMENSIONAL ASSIGNMENT PROBLEMS 

III.l. Introduction 

The ever- increasing demand for sensor surveillance 

systems is to accurate track and identify objects in real- 

10 time, even for dense target scenarios and in regions of high 
track contention. The use of multiple sensors, through more 
varied information, has the potential to greatly improve 
target state estimation and identification. However, to take 
full advantage of the available data, a multiple frame data 

15 association approach is needed. The most popular such 
approach is an enumerative technique called multiple 
hypothesis tracking (MHT) . As an enumerative technique, MHT 
can be too computationally intensive for real-time needs 
because this problem is NP-hard. A promising approach is to 

20 utilize the recently developed Lagrangian relaxation 
algorithms (K. P . Pattipati , S. Deb, and Y . Bar-Shalom. A 
multisensor-multitarget data association algorithm for 
heterogeneous sensors. IEEE Transactions on Aerospace and 
Electronic Systems, 29, No. 2:560-568, April 1993; A.B. Poore 

25 and N.Rijavec. Partitioning multiple data sets: 
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multidimensional assignments and lagrangian relaxation, in 
quadratic assignment and related problems, P.M. Pardalos and 
H.Wolkowicz, editors, DIMACS series in Discrete Mathematics 
and Theoretical Computer Science, 16:25-37, 1994; A.J. 
5 Robertson III. A class of lagrangian relaxation algorithms for 
the multidimensional problem. Ph.D. Thesis, Colorado State 
University, Ft. Collins, CO, 1995 J for the multidimensional 
assignment problem; however, there are many other potentially 
good approaches to these assignment problems such as LP 

10 relaxation combined with an interior point method, GRASP, and 
parallel izat ion. 

These data association problems are fundamentally 
important in tracking. Thus, our first objective is to 
present an overview of the tracking problem and the context 

15 within which these data association problems fit. Following 
a brief review of the probabilistic framework for data 
association, the second objective is to formulate two models 
for track initiation and maintenance as multidimensional 
assignment problems. This formulation is based on a window 

20 moving over the frames of data. The first and simpler model 
uses the same window length for track initiation and 
maintenance, while the second model uses different lengths. 
As the window moves over the frames of data, one creates a 
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sequence of these multidimensional assignment problems and 
there is an overlap region of frames common to adjacent 
windows. Thus, given the solution of one problem, one can 
develop a warm start for the solution to the next problem in 
5 the sequence, as shown hereinafter. Such information is 
critically important to the design of real-time algorithms. 

III. 2. Overview of the Tracking Problem 
Tracking and data fusion are complex subjects of 

10 specialization that can only be briefly summarized as they are 
related to the subject of this paper. The processing of track 
multiple targets using data from one or more sensors is 
typically partitioned into two stages or major functional 
blocks, namely, signal processing and data processing (Y.Bar- 

15 Shalom. Multitarget-Multisensor Tracking: Advanced 
Applications. Artech House, Dedham, MA, 1990; Y. Bar-Shalom and 
T.E. Fortmann. Tracking and Data Association. Academic Press, 
Boston, MA, 1988; S.S. Blackman. Multiple Target Tracking with 
Radar Applications. Artech House, Dedham, MA, 1986) . The 

2 0 first stage is the signal processing that takes raw sensor 
data and outputs reports. Reports are sometimes called 
observations, threshold exceedances, plots, hits, or returns, 
depending on the type of sensor. The true source of each 
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report is usually unknown and can be due to a target of 
interest, a false signal, or persistent background objects 
that can be moving in the field of view of the sensor. 

Two principal functions of data processing are tracking 
5 and discrimination or target identification. The 
discrimination or target identification function distinguishes 
between tracks that are for targets of interest and other 
tracks such as those due to background. Also, if there is 
enough information, each track is classified as to the type of 

10 target it appears to represent. Target identification and 
discrimination will not be discussed further in this paper 
except to comment here on the use of attributes (including 
identification information) in the data association. 

There are two types of attributes or features, namely, 

15 discrete valued variables and continuous valued variables. 
Available attributes of either or both types should be used in 
computing the log likelihoods or scores for data association. 
In discussing tracking in this paper, the attributes will not 
be mentioned explicitly but it should be understood that some 

2 0 reports and some tracks may include attributes information 
that should be and can be used in the track processing (both 
data association and filtering) if it is useful. 
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III. 2.1. Tracking Functions 
The tracking function accepts reports for each frame of 
data and constructs or maintains a target state estimate, the 
variance-covariance matrix of the state estimation error, and 
5 refined estimates (or probabilities) of the target attributes. 
The state estimate typically includes estimates of target 
position and velocity in three dimensions and possibly also 
acceleration and other variables of interest. 

The tracking function is typically viewed in two stages, 
10 namely, data association and filtering. Also, the process of 
constructing new tracks, called track initiation, is different 
from the process of updating existing tracks, called track 
maintenance . 

In track maintenance, the data association function 
15 decides how to relate the reports from the current frame of 
data to the previously computed tracks. In one approach, at 
most one report is assigned to each track, and in other 
approaches, weights are assigned to the pairings of reports to 
a track. After the data association, the filter function 
2 0 updates each target state estimate using the one or more (with 
weights) reports that were determined by the data association 
function. A filter commonly used for multiple target tracking 
and data fusion is the well known Kalman filter (or the 
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extended Kalman filter in the case of nonlinear problems) or 
a simplified version of it (Y. Bar-Shalom. Multitarget- 
Multisensor Tracking: Advanced Applications . Artech House, MA, 
1990; Y. Bar-Shalom and T.E. Fortmann. Tracking and Data 
5 Association. Academic Press, Boston, MA, 1988; S.S. Blackman. 
Multiple Target Tracking with Radar Applications . Artech 
House, Dedham, MA, 1986) . 

In track initiation, typically a sequence of reports is 
selected (one from each of a few frames of data) to construct 

10 a new track. In track initiation, the filtering function 
constructs a target state estimate and related information 
based on the selected sequence of reports. The new track is 
later updated by the track maintenance processing of the 
subsequent frames of reports. 

15 In some trackers, there is also a tentative tracking 

function for processing recently established tracks until 
there is enough confidence to include them in track 
maintenance. While for simplicity of presentation, tentative 
tracking will not be included in the processing discussed in 

2 0 this paper, the techniques discussed could readily include one 
or more tentative tracking functions. 

There have been numerous approaches developed to perform 
the data association function. Since optimal data association 
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is far too complex to implement, good but practical sub- 
optimal approaches are pursued. Data association approaches 
can be classified in a number of ways. One way to classify 
data association approaches is based on the number of data 
5 frames used in the association process (O.E. Drummond. 
Multiple sensor tracking with multiple frame, probabilistic 
data association. In Signal and Data Processing of Small 
Targets, SPIE Proceedings, volume 2561, pages 322-336, 1995) . 
In single frame data association for track maintenance, 

10 "hard decisions" are made on the assignment of reports to the 
tracks. Some single frame approaches include: individual 
nearest neighbor, PDA, JPDA, and global nearest neighbor 
(sequential most probable hypothesis tracking) , which uses a 
two-dimensional assignment algorithm (Y. Bar- Shalom. 

15 Multitarget-Multi sensor Tracking: Advanced Applications . 
Artech House, MA, 1990; Y. Bar-Shalom and T.E. Fortmann. 
Tracking and Data Association. Academic Press, Boston, MA, 
1988; S.S. Blackman. Multiple Target Tracking with Radar 
Applications. Artech House, Dedham, MA, 1986) . In most single 

2 0 frame data association approaches, only one track per object 
is carried forward to be used in processing the next frame of 
reports . 
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Multiple frame data association is more complex and 
frequently involves purposely carrying forward more than one 
track per target to be used in processing the next frame of 
reports. By retaining more than one track per target, the 
5 tracking performance is improved at the cost of increased 
processing load. The best known multiple frame data 
association approach is an enumerative technique called 
multiple hypothesis tracking (MHT) which enumerates all the 
global hypothesis with various pruning rules. 

10 As shown hereinafter, the log likelihood of the 

probability of a global hypothesis can be decomposed into the 
sum of scores of each track that is contained in the 
hypothesis. As a consequence, the most probable hypothesis 
can be identified with the aid of a non- enumerative algorithm 

15 for the assignment so formulated. 

III. 2. 2. Interface Between Track Initiation 
And Track Maintenance 

There are a number of ways that track initiation and 

track maintenance functions can interact (O.E. Drummond. 

2 0 Multiple target tracking with multiple frame, probabilistic 

data association. In Signal and Data Proceedings, volume 1954, 

pages 394-408, 1993) . Two ways that are especially pertinent 

to the methods of this paper will be discussed in this 

section. 
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One approach is to treat these two functions 
sequentially, by first assigning reports to tracks and then 
using the information that remains to form new tracks. A 
better approach is to integrate both processes where 
5 initiating tracks compete for the new frame of reports on an 
equal basis, i.e., at the same time. 

In the integrated approach, the data association for 
track initiation and track maintenance processing are combined 
and conducted simultaneously. One assignment array is created 

10 that includes the track scores for potential tracks for both 
track initiation and track maintenance. The first dimension 
of the assignment problem includes only all the tracks either 
created or updated in the processing of the prior frames of 
reports. Each of the remaining dimensions accommodates one 

15 frame of reports. 

After the assignment algorithm finds a solution, each of 
the previously established tracks are updated with the report 
assigned to it from the second dimension of the assignment 
array. The remaining reports in the second dimension that 

2 0 were in the assignment solution are firmly established as 
track initiators. The unassigned reports in the second 
dimension of the assignment array are discarded as false 
signals. The processing of the next frame of reports repeats 
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this process using these updated tracks and the newly 
identified initiators in the first dimension of the cost array 
for processing the new frame of reports. Details of this 
approach are discussed hereinafter. 
5 The integrated approach just discussed, uses the same 

number of frames of reports for both track initiation and 
track maintenance. The goal in using the mult i -dimensional 
assignment algorithm is to provide improved performance while 
minimizing the amount of processing required. Typically, 

10 track initiation will benefit from more frames of reports than 
will track maintenance. Thus, a second approach that 
integrates track initiation and track maintenance is discussed 
hereinafter, wherein the number of frames of reports is not 
the same for these two functions. This is a novel approach 

15 that is introduced in this paper. 

III. 2. 3. Multiple Sensor Processing 
There are many advantages and many ways of combining data 
from multiple sensors. There are also many ways of 
2 0 categorizing the different algorithmic architectures for 
processing data from multiple sensors. One approach outlines 
four generic algorithmic architectures (O.E. Drummond. 
Multiple sensor tracking with multiple frame, probabilistic 
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data association. In Signal and Data Processing of Small 
Targets, SPIE Proceedings , volume 2561, pages 322-336, 1995) . 
Two of these generic architectures are especially pertinent to 
this paper and are summarized briefly. 
5 In the Centralized Fusion algorithmic architecture, 

reports are combined from the various sensors to form global 
tracks. This algorithmic architecture is also called Central 
Level Tracking, Centralized Algorithmic Architecture, or 
simply the Type IV algorithmic architecture. In track 

10 maintenance, for example, if a single frame data association 
is used, then for data association a frame of reports from one 
sensor is processed with the latest global tracks; then the 
global tracks are updated by the filter function. After the 
processing of this frame of reports is completed, a frame of 

15 the reports from another (or the same) sensor is processed 
with these updated global tracks. This process is continued 
as new frames of data become available to the system as a 
whole . 

In Centralized Fusion, using the multi-dimensional 
2 0 assignment algorithm for the data association with multiple 
sensors is similar to the processing of data from a single 
sensor. Instead of using multiple frames of reports from a 
single sensor, the multiple frames of reports come from 
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multiple sensors. The frames of reports from all the sensors 
are ordered based on the nominal time each frame of reports is 
obtained and independent of which sensor provided the reports. 
In this way the approaches discussed in this paper can be 
5 extended to process reports from multiple sensors using the 
Centralized Fusion algorithmic architecture. 

The second pertinent algorithmic architecture is Track 
Fusion. This approach is also called the Hierarchical or 
Federated algorithmic architecture, sensor level tracking or 
10 simply the Type II algorithmic architecture. In track 
maintenance, for example, a processor for each sensor computes 
single sensor tracks; these tracks are then forwarded to the 
global tracker to compute global tracks based on data from all 
sensors . 

15 After the first time a sensor processor forwards tracks 

to the global tracker, then subsequent tracks for the same 
targets are cross-correlated with the existing global tracks. 
This track-to-track cross-correlation is due to the common 
history of the current sensor tracks and the tracks from the 

20 same sensor that were forwarded earlier to the global tracker. 
The processing must take this cross-correlation into account 
and there are a number of ways of compensating for this- cross - 
correlations. One method for dealing with this cross- 
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correlation is to decorrelate the sensor tracks that are sent 
to the global tracker. There are a variety of ways to achieve 
this decorrelation and some are summarized in recent paper 
(O.E. Drummond. Feedback in track fusion without process 
5 noise. In Signal and Data Processing of Small Targets, SPIE 
Proceedings, volume 2561, pages 369-383, 1995) . 

Once the sensor tracks are decorrelated they can be 
processed by the global tracker in almost the same way as 
reports are processed. In the case of track fusion the 

10 association process is referred to as track (or track-to- 
track) association rather than data (or report-to-track) 
association. If the sensor tracks are decorrelated, the 
global tracker can process the tracks from the various sensor 
processors in much the same way that the global tracker of 

15 Centralized Fusion processes reports. Accordingly the methods 
described in this paper can be readily extended to processing 
data from multiple sensors using either the Centralized Fusion 
or Track Fusion or even a hybrid combination of both 
algorithmic architectures . 

20 

III. 3. Formulation of the Data Association Problem 
The goal of this section is to briefly outline the 
probabilistic framework for the data association problems 
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presented in this work. The technical details are presented 
elsewhere (A.B. Poore. Multidimensional assignment formulation 
of data association problems arising from multitarget tracking 
and multisensor data fusion. Computational Optimization and 
5 ApplicationSjr 3:27-57, 1994). The data association problems 
for multisensor and multitarget tracking considered in this 
work are generally posed (A.B. Poore. Multidimensional 
assignment formulation of data association problems arising 
from multitarget tracking and multisensor data fusion. 
10 Computational Optimization and Applications, 2:21-51, 1994) as 
that of maximizing the posterior probability of the 
surveillance region (given the data) according to 

Maximize {P{F=y\ Z ml ) \ y e r*}, (3.1) 
where Z N+1 represents N+1 data sets, y is a partition of indices 

15 of the data (and thus induces a partition of the data) , r* is 
the finite collection of all such partitions, r is a discrete 
random element defined on r*, and P(r=Yl zN+1 ) is the posterior 
probability of a partition y being true given the data Z M+1 . 
The term partition is defined below; however, this framework 

2 0 is currently sufficiently general to cover set packings and 
coverings . 
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Consider N+l data sets Z(k) (k=l, . . . , N+l) , each consisting 
of M k actual reports and a dummy report , and let denote the 
cumulative data set defined by 

5 Z(k) and Z ml = {Z(l) ,...,Z(N+1)}, (3 . 2 ) 

respectively. (The dummy report z 0 serves several purposes in 
the representation of missing data, false reports, initiating 

10 tracks, and terminating tracks (A.B. Poore . Multidimensional 
assignment formulation of data association problems arising 
from multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications, 3:27-57, 1994).) 
In multisensor data fusion and multitarget tracking the data 

15 sets Z(k) may represent different classes of objects, and each 
data set can arise from different sensors. For track 
initiation the objects are reports that must be partitioned 
into tracks and false alarms. In our formulation of track 
maintenance, which uses a moving window, one data set will be 

2 0 tracks and remaining data sets will be reports which are 
assigned to existing tracks, as false reports, or to 
initiating tracks. 

We specialize the problem to the case of set partitioning 
(A.B. Poore. Multidimensional assignment formulation of data 

25 association problems arising from multitarget tracking and 
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multisensor data fusion. Computational Optimization and 
Applications, 3:27-57, 1994) defined in the following way. 
Define a "track of data" as iz, 1 , zf +1 \ where each i. can 
assume zero or nonzero values. A partition of the data will 
refer to a collection of tracks of data wherein each report 
occurs exactly once in one of the tracks of data and such that 
all data is used up; the occurrence of dummy report is 
unrestricted. The reference partition y° is that in which all 
reports are declared to be false. 

Next, under appropriate independence assumptions one can 
show (A.B. Poore. Multidimensional assignment formulation of 
data association problems arising from multitarget tracking 
and multisensor data fusion. Computational Optimization and 
Applications, 3:27-57 , 1994) 

?(r= Y \Z^) n 

P (r = Y°| Z N+1 ) e Y W ^ * ; 

L i x .i N+1 is a likelihood ratio containing probabilities for 
detection, maneuvers, and termination as well as probability 
density functions for report errors, track initiation and 
termination. Define 



c = -InL. 

-Z-i -^M-ul -2-1 



/ 1, if(z. .. . ) are assigned to as a track, 

z ~ \ *i (3 4) 

2 i z n+i l o, otherwise. 
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Then, recognizing that - In 



P(r = y|z N ) 
p(r = y°|2 n+1 ) J 



the problem (3.1) can be expressed as the (N+l) -dimensional 
assignment problem 



Minimize • ■ • 2^, °- ■ z - ■ 

i =0 1=0 

Subject To E... E z i =1, 1=1,...,**,, 



for i k = l,...,M k and k = 2,...,N, 



0.5; 



5 where c 0 0 is arbitrarily defined to be zero. Here, each group 
of sums in the constraints represents the fact that each non- 
dummy report occurs exactly once in a "track of data." One 
can modify this formulation to include mult i -assignments of 
one, some, or all the actual reports. The assignment problem 
10 (3.5) is changed accordingly. For example, if z* is to be 
assigned no more than, exactly, or no less than nf times, then 
the "=1" in the constraint (3.5) is changed to "<, =, > nf , " 
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respectively. In making these changes, one must pay careful 
attention to the independence assumptions, which need not be 
valid in many applications. Expressions for the likelihood 
ratios can be found in the work of (A.B. Poore . 
5 Multidimensional assignments and multitarget tracking, 
partitioning data sets, I.J. Cox, P. Hansen, and B. Julesz, 
editors. DIMACS Series in Discrete Mathematics and Theoretical 
Computer Science, American Mathematical Society, Providence, 
JR. I., 19:169-198, 1995) and the references therein. 

10 

III. 4. Track Initiation and Maintenance 
In this section we explain two multiframe assignment 
formulations to the track initiation and maintenance problem. 
The continued use of all prior information is computationally 

15 intensive for tracking, so that a window sliding over the 
frames of reports is used as the framework for track 
maintenance and track initiation within the window. The 
objectives are to describe an improved version of a simple 
method and then to put this into a more general framework in 

20 which track initiation and maintenance have different length 
moving windows . 
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III. 4.1. The First Approach to Track 
Maintenance and Initiation 

The first method as explained in this section is an 

improved version of our first track maintenance scheme (A.B. 

5 Poore. Multidimensional assignment formulation of data 

association problems arising from multi target tracking and 

multisensor data fusion. Computational Optimization and 

Applications, 3:27-57, 1994) and uses the same window length 

for track initiation and maintenance after the initialization 

10 step. The process is to start with a window of length I\T+1 
anchored at frame one. In the first step there is only one 
track initiation in that we assume no prior existing tracks. 
In the second and all subsequent frames, there is a window of 
length N anchored at frame k plus a collection of tracks up to 

15 frame k. This window is denoted by [k, Jc+1,..., k+N} . The 
following explanation of the steps is much like mathematical 
induction in that we explain the first step and then step k to 
step Jc+1. 

Track Maintenance and Initiation: Step 1. Let 

be an enumeration of all those zero-one variables in the 
solution of the assignment problem (3.5) (i.e., z .. . =1) 

1 l 1 2' 1 M+1 

excluding all the false reports in the solution (i.e., all 
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those zero-one variables with exactly one nonzero index) and 
zero-one variables in the solution for which {i 1 ,i 2 ) = (0,0). 
(The latter can correspond to tracks that initiate on frames 
three and higher.) These denote our initial tracks. 
5 Consider only the first two index sets in this 

enumeration (3.6) and add the zero index 1 2 =0 with the 
corresponding values of ± 1 and i 2 being zero. Thus, the 
enumeration is now {i 1 (1 2 ) , i 2 (1 2 ) . The notation 

T 2 {1 2 ) = (^{^j ' z i <i ) ) will be used for this pairing. Suppose 
10 now that the next data set, i.e., the (i\T+2) th set, is added to 
the problem. 

To explain the costs for the new problem, one starts with 
the hypothesis that a partition yET* being true is now 
conditioned on the truth of the pairings on the first two 
15 frames being correct. Corresponding to the sequence 
{ T 2 ^2)' z i 3 ' z i N+2 \ ' the likelihood function is then given by 

L= II L , where 

{r 2 u 2 >, v ...,* iiM }6Y ( 3. 7) 

Li = T T 

L T 2 (0) =1 ' 

Next, define the cost and the corresponding zero-one variable 

20 
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C L = -lni, 



z . . 

J -2 X 3"* 1 W+2 



1, if j^y is assigned to Track T 2 {1 2 ) , 
0, otherwise, 



(3.8) 



respectively. Then the track maintenance problem Maximize 
{lylyEr*} can be formulated as the following multidimensional 
assignment 



Minimize Jj ^ . .... . ... . 

1 2 =0 i 3 = 0 i N+2 = 0 2 3 w+2 2 3 w+2 



Subject to: ^ (3.9) 

Z-^ .2^ ^7 i „j = 1/ 1 9 -1 F ...,L r 



V 0 W° 

L, M, 



2, 2-, • • • Z li...i = 1 r i 3 = l r -,M r 

2-i 2—i ••• 2—i '•' 2—i z i ± ... 7 ~l/ 

1 2 =0 i 3 =0 Vi=° W° 23 N * 2 

for i =1,...,M and p=4,... f tf+2-1, 



z ii...i ^{0,1} for all 1 



Track Maintenance and Initiation: Step k. At the 

beginning of the Jc* 1 step, we solve the following (Jff+i) - 
dimensional assignment problem. 
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Minimize £ • • • £ _ z ± ...i (3.10) 



Subject to: 



i =n i -n J -k ± k+i — ± k+N K K 

1 k+l u k+N 

X, 2-/ ••• £ ^ i =1/ i^ +1 = l,-..M , 

2s 2s * • ■ 2s 2s * • • 2s z i i i " ^ 

i=n i =n i =n-i =n i =n * £ + i"* k+N 
for i =1,.../M_ , andp=Jc+2 f ... / N+ic-l A 



2. 2s .... 2. n ^ w -i^ =1 ' W 1 -.,^ 

z 7 .... eiO,lf for all 1. , i. +1 , 



where for the first step I x and L x are replaced by i 2 and M x , 
respectively. The first index I k in the subscripts corresponds 
to the sequence of tracks {r.(IJP , where 

\ K K J 1^ = 0 

T k^ 1 k ) = { z i"i Z ± k ^ 1 k i } is a track of data from the solution 

10 of the problem prior to the formulation of (3.10). If k=l, 
then it is just the first data set in frame one. 

Basic Assumption. Suppose problem (3.10) has been solved 
and let the solution, i.e., those zero-one variables equal to 
one, be enumerated by 

with the following exceptions. 
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a. All zero-one variables for which [l k , i k+1 ) - (0, 0) are 
excluded. Thus, tracks that initiate on frames after the 
(k+l) th are not included in the list. 

b. All zero-one variables whose subscripts have the form I k =0 
and exactly one nonzero index in the remaining indices 
{ i k+1 r i k+N ) are excluded. These correspond to false 
reports on frames p = k+l,...,k+N. 

c. All variables z for which (i k+1 ,... r i k+N ) = (0, 0,..., 0) 
and i Jt (I jt )=0 are excluded from (1.11). In other words, 
the reports on the last N+l frames in the 

{ T k ^k) ' z K+i f * * " ' Z K^ are dummy. An Y solution with 

this feature corresponds to a terminated track. 
Given the enumeration (3.11), one now fixes the assignments 
on the first two index sets in the list (3.11). The zero 
index l k+1 = Q is added to the enumeration to specify 
(^(0) , i k+1 (Q))= (0, 0) that is used to represent false reports and 
tracks that initiate on frame k+2 or later, so that the 
enumeration (3.11) is now fe(I,J,i, +1 (I w ))}j;; B0 . 

Then, for l k+1 = l , L k+1 , the 1^ such track is denoted by 

^1(^1)= KtV 1 *^)' z iil {1 ^i ] } and the (N+l) -tuple 
{ r *+i( 2 * + i> < z Zl<-' z Z^\ wil1 denote a track T k+1 (l k+1 ) plus a set 
of reports \ z i k+2 ''~, Z i k+1+N }' actual or dummy, that are feasible 
with the track 2^(1^) . The (N+l) -tuple {^{0), z^,„.,z^} 
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will denote a track that initiates in the sliding window, 
i.e., on subsequent frames. A false report in the sliding 
window is one with all but one non-zero index i p for some 
p=k+2,...,k+l+N in the -tuple {^(0), z Zl' ' z Z!x + u\ ' 

5 The corresponding hypothesis about a partition y E r* 

being true is now conditioned on the truth of the L k+1 tracks 
existing at the beginning of the N-frame window. (Thus the 
assignments prior to this sliding window are fixed.) The 
likelihood function is given by 

L = IT L , where 



10 



, L k+2 



1. 



(3.12; 



Next, define the cost and the zero-one variable by 

{ / * 2 k 1 1 (3.13) 

z . . if \ z Zl> " ' Z W J is assigned to , 

-WW-Wm 0, otherwise, 



15 



respectively, so that the track extension problem, which was 
originally formulated as Maximize {L y \yGT*j f can be expressed 
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as exactly the same multi-dimensional assignment in (3.4) with 
k replaced by k+1 . Thus, we do not repeat it here. 



III. 4. 2. The Second Approach to Track 
Maintenance and Initiation 



In the first approach the window lengths were the same 



for both track maintenance and initiation. 



This can be 



inefficient in that one can usually use a shorter window for 
track maintenance than for track initiation. This section 

10 addresses such a formulation. 

The General Step k. To formulate a problem for track 
initiation and maintenance we consider a moving window 
centered as frame k of length T-k7+1 denoted by 
[&-!,..., k,...,k+J] . In this representation, the window length 

15 for track maintenance is J and that for initiation, I+J+l. The 
objective will be to explain the situation at this center and 
then the move to the same length window at center k+1 denoted 
by [k+1 -I,..., Jc+1,..., k+l+J] , i.e., by moving the window to the 
right one frame at time. The explanation from the first step 

20 follows hereinafter. 



The notation for a track of data is 




(3 .14) 
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where the index l k is used for an enumeration of those reports 



the sequence of reports belonging to track T k (l k ) but 
restricted to frames prior to and including p. Thus, 



Given this notation for the tracks and partition of the 
data in the frames { k-I, k, k+ J} , L t k a k ) will denote the 
accumulated likelihood ratio up to and including frame p(p<k) 
for a track that is declared as existing on frame Jc as a 
10 solution of the assignment problem. In this notation, the 
likelihood for T Pfk i ^ 1 k ) and that of the association of 



paired together. We also use the notation T pfk {l k ) to denote 



5 




for any p<k-l. 



(3.15) 




with track T k (l k ) is given by 




for any p < k- 1, 



(3.16) 



15 respectively. 
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z 



The cost for the assignment of j z i^ * r z ±^ j to trac k 
T k (l k ) and the corresponding zero-one variable are given by 



c 7 7 7 - = ™ln [L_ fl ,L. . ] 

= " ln[ \^ i ^ i vi( I P-" i ^vw-iJ 

(3.17) 

t ) 1, if" |z*^,... f z^*^ | is assigned to track T k (l k ) , 

0, otherwise, 



respectively. Likewise, for costs associated with the false 
reports on the frames k-I to k and as associated with 
10 jz^,...,^^! and the corresponding zero-one variables are 
given by: 



i = -ln[i. 



k-l ^k^k+1 

( 1, if {z, z. 



(3 .18) 



15 



are assigned as a track, 
0, otherwise. 



For the window of frames in the range {k-I,...,k} , the easiest 
explanation of the partition of the reports is based on the 
definitions 
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p,k 



' i \zf belongs to one of the tracks 

listed on frame k, i.e., to one \r (3.19) 
of the TA1,) for some I =1,...,I. 



so that all of those reports not used in the tracks T k ( 1 k ) on 
frame p are put into the set of available reports F Pfk for the 
formation of tracks over the entire window. Thus, we 
formulate the assignment problem as 



Minimize £ £ ^ c[ i i i z[ ± t ± 

ptf-; tfi ptk r trt +1 t^Q ^-r^JrW^+jr ^r-V k+ v^ k 

+ z2 X, 22 °i i ± z l ± 



1^=1 r=Jt + l i r =0 ^i-^J 



Subject To £ E ? £^i F ii i =1 (3-20) 

{p=*-I,p*qr} i €F ptk r=k+l i r =0 *~ r * * +1 w 

for i g 6F p/Jt \{0} a/2d q=k-I,».,k, 

p=*-I i p €F pfk {r=k+l,r*q} i q =0 k ~ r k * +1 * +J 

for i g =l,.„, Af and q= k+1 , ... , k+J, 

Z ±r-Vn-^ €{0 ' 1} f0raU Vl-V^l^-^ 



r=£+l i=0 



Z V, + l-A +J =1 ' i q =1 '-' M q 
l k =l {r=k*l,r*q} i r =0 * * +1 k * J q 9 

for q=k+l, ... r k+J, 

\^ k J {0 ' 1} f ° r 511 Vl'W-Vl- 



140 



CSUR.01USRI 

Basic Assumptions . Suppose problem (1.2 0) has been solved 
and let the solution, i.e. those zero-one variables equal to 
one, be enumerated by 

T 

(3.21) 



f F 



with the following exceptions. 

a. All zero-one variables in the second list for which 

( 1 k - I f-f i k r i k+1 ) ~ (0 f 0, 0) are excluded. Thus, tracks that 
initiate on frames after the (Jc+1) th are not included in 
the list. 

b. All false reports are excluded, i.e., all zero-one 
variables in the second list whose subscripts have 
exactly one nonzero index. 

c. All variables z i k i k+lm i k+lf for which ( i k+1 , - r i k+N ) = ( 0 , 0 , ... , 0 ) 
and Z{p) nT k {l k ) =0 f or p=k-K,...,k where K >0 is user 
specified. Thus the track T k (l k ) ± s terminated if it is 
not observed over K+J+l frames. 

Given the enumeration (3.20) , one now fixes the assignments on 
the all index sets up to and including the (Jc+1) th index sets. 
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for l k+1 = l,...,L k ^, 



T (1 ) =<, 



) 



(3.22) 



for I j t+ i = £ * + i + 1 ^--/ i Jk+ i- 




Thus one can now formulate the assignment problem for the 
next problem exactly as in (3.20) but with k replaced by k+1. 
Thus, we do not repeat it here. 

The Initial Step. Here is one explanation for the 
5 initial step. First, assume that N=X+J. In this case, we 
start the track initiation with a solution of (3.5). Let 



be an enumeration of the solution set of (3.5), i.e., those 
zero-one variables z i 1 (i 2 )± 2 (i 2 )...i N+1 a 2 ) =1 , including z 00 0 =1 

10 corresponding to 1 2 =0, but excluding all those zero-one 
variables that are assigned to one and correspond to false 
reports (i.e., there is exactly one nonzero index in the 
subscript of z ± 1 i 2 -.i N+1 ) , all those zero-one variables that are 
assigned to one and correspond to tracks that initiate on 

15 frames higher than 1+2. Then we fix the data association 
decisions corresponding to the reports in our list of tracks 
prior to and including frame k+l=I+2. This defines the k for 




± 2 
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problem (3.5) and one can then continue the development by 
adding a frame to the window as in the general case . 

If T+J" > N, then one possibility is to start the process 
with N+l frames, and assuming J<N f proceed as before replacing 

5 I by N-J for the moment, and continue to add frames without 
lopping off the first frame in the window until reaches a 
window of length J+J+l. Then we proceed as in the previous 
paragraph . 

If J+J < N, then one can solve the track initiation 
10 problem (3.5), formulate the problem with the center of the 
window at k+1 = N+l- J, enumerate the solutions as above, and 
lop off the first N-J -I frames. Then, we proceed just as in 
the case I+J=N. 

III. 5. CONCLUDING COMMENTS 
15 A primary objective in this work has been to demonstrate 

how multidimensional assignment problems arise in the tracking 
environment. The problem of track initiation and maintenance 
has been formulated within the framework of a moving window 
over the frames of data. The solution of these NP-hard, 
2 0 noisy, large scale, and sparse problems to the noise level in 
the problem is fundamental to superior track estimation and 
identification. Thus, one must utilize the special structure 
in the problems as well take advantage of special information 
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that is available. Since these moving windows are 

overlapping, there are some algorithm efficiencies that can be 
identified and that take advantages of the overlap in the 
windows from one frame of reports to the next. Here is an 
example of the use of a primal solution of one problem to warm 
start the solution of the next problem in the sequence. 

Suppose we have solved problem (3.10) and have enumerated 
all those zero-one variables in the solution of (3.10) as in 
(3.11). Add the zero index l k+1 =Q , so that the enumeration is 

{( h < h+i > ri k+1 ( I w ),-.., i k+N d k+1 ) )}^ +1 _ 0 . (3.23) 

With this enumeration one can define the cost by 

o . ~ c (324) 

and the two-dimensional assignment problem 

® 2 = Minimize £J £ c 2 2 ± z\ i =VAz 2 ) 

* 1 -C\ -t -n ^A+l^A+l+w ± k+i 1 k+l+N 2 

X *+1~ U 2 Jt+l+Af U 



Subject to X < lWw =l, l k+1 = l,...,L k+1 

•t-k+l+N'V 

2 

-^4.1 ~~ ^ 



(3.25) 



■Jk+i" 

z L^j^>^ for * u -w 



Let w be an optimal or feasible solution to this two- 
dimensional assignment problem and define 
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z . . - 



1, if [i k+1 ,...,i k+N ) =(i k + 1 (l k+1 ),...,i k+N (l k+1 )) 

and = 1 f ° r SOme Wl'-'-^i {3 26) 

or if (l w ,i wtl( ] = (0,0) 

0, otherwise. 



This need not satisfy the constraints in that there are 
5 usually many objects left unassigned. Thus, one can complete 

the assignment by using the zero-one variables in (3.10) with 

k replaced by Jc+1 with exactly one nonzero index corresponding 

to any unassigned object or data report. 

For the dual solutions, the multipliers arising from the 
10 solution of the two-dimensional assignment problem (3.25) 

corresponding to the second variable, i.e., | U * +1+W P 1+S m 

These are good initial values to use in a relaxation scheme 
(A.B. Poore and N.Rijavec. Partitioning multiple data sets: 

multidimensional assignments and Lagrangian relaxation, in 
15 quadratic assignment and related problems; P.M. Pardalos and 

H. Wolkowicz, editors, DIMACS series in Discrete Mathematics 

and Theoretical Computer Science, 16:25-37, 1994; A.J. 

Robertson III. A class of lagrangian relaxation algorithms for 

the multidimensional assignment problem. Ph.D. Thesis, 
20 Colorado State University, Fr. Collins, CO, 1995) . Finally, 

note that one can also develop a warm start for problem (3.20) 

in a similar fashion. 
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IV. REVIEW OF NEW RELAXATION SCHEMES 
IV. 1. Introduction 
Multidimensional assignment problems govern the central 
problem of data association in multisensor and multitarget 
tracking, i.e., the problem of partitioning observations from 
multiple scans of the surveillance region and from either 
single or multiple sensors into tracks and false alarms. This 
fundamental problem can be stated as (A.B. Poore . 
Multidimensional assignments and multitarget tracking, 
partitioning data sets, I.J. Cox, P. Hansen, and B. Julesz, 
editors, DIMACS Series in Discrete Mathematics and Theoretical 
Computer Science, American Mathematical Society, Providence, 
R.I., 19:169-198, 1995; A.B. Poore. Multidimensional 
assignment formulation of data association problems arising 
from multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications, 2:21-51, 1994) 
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... £ c, . z. . 

i i=0 i M=0 i' » 

Subject to: (4 .D 

£ • • • £ z =1, i=l f ...,M, 

i 2 =0 i N =0 ^ ^ 1 

• • • 2!^ z^ • • • Z^ z i ...± =1 

J i =0 Vr° Vi=° V° 1 " 
for i p =l,...,M p a^d p=2, A7-1, 

z ^ =1 ' i » =1 "- M »' 

where c 0 ... 0 is arbitrarily defined to be zero and is included 
for notational convenience. One can modify this formulation 
to include multi- assignments of one, some, or all of the 
actual reports. The zero index is used in representing 
missing data, false alarms, initiating and terminating tracks 
as described in (A.B. Poore. Multidimensional assignments and 
multitarget tracking, partitioning data sets, I.J. Cox, P. 
Hansen, and B. Julesz, editors, DIMACS Series in Discrete 
Mathematics and Theoretical Computer Science, American 
Mathematical Society, Providence, R.I., 19:169-198, 1995; A.B. 
Poore. Multidimensional assignment formulation of data 
association problems arising from multitarget tracking and 
multisensor data fusion. Computational Optimization and 
Applications, 3:27-57, 1994). In these problems, we assume 
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that all zero-one variables z . . with precisely one nonzero 

J -1'"~ L N 

index are free to be assigned and that the corresponding cost 
coefficients are well-defined. (This is a valid assumption in 
the tracking environment (A.B. Poore. Multidimensional 
assignments and multitarget tracking, partitioning data sets, 
I.J. Cox, P. Hansen, and B. Julesz, editors, DIMACS Series in 
Discrete Mathematics and Theoretical Computer Science, 
American Mathematical Society, Providence, R.I., 19:169-198, 
1995; A.B. Poore. Multidimensional assignment formulation of 
data association problems arising from multitarget tracking 
and multi sensor data fusion. Computational Optimization and 
Applications, 3:27-57, 1994).) Although not required, these 
cost coefficients with exactly one nonzero index can be 
translated to zero by cost shifting (A.B. Poore and N. 
Rijavec. A lagrangian relaxation algorithm for 
multidimensional assignment problems arising form multitarget 
tracking. SIAM Journal of Optimization, 3, No. 3:544-563, 
1993) without changing the optimal assignment. Finally, this 
formulation is of sufficient generality to include the 
symmetric problem and the asymmetric inequality problem (A.J. 
Robertson III. A class of lagrangian algorithms for the 
multidimensional assignment problem. Ph.D. Thesis, Colorado 
State University, Fr. Collins, CO, 1995) . 
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The data association problems in tracking that are 
formulated as Equation (4.1) have several characteristics. 
They are normally sparse, the cost coefficients are noisy and 
the problem is NP-hard (M.R. Galey and D.S. Johnson. Computers 
and Intractability. W.H. Freeman and Company, San Francisco, 
CA, 1979), but it must be solved in real-time. The only known 
methods for solving this NP-hard problem optimally are 
enumerative in nature, with branch-and-bound being the most 
efficient; however, such methods are much too slow for real- 
time applications. Thus one must resort to suboptimal 
approaches. Ideally, any such algorithm should solve the 
problem to within the noise level, assuming, of course, that 
one can measure this noise level in the physical problem and 
that the mathematical method provides a way to decide if the 
criterion has been met. 



There are many algorithms that can be used to construct 
suboptimal solutions to NP-hard combinatorial optimization 
problems. These include greedy (and its many variants), 
relaxation, simulated annealing, tabu search, genetic 
algorithms, and neural network algorithms (C.H. Papadimitriou 
and K.Steiglitz. Combinatorial Optimization: Algorithm and 
Complexity. Prentice-Hall, Inc., Engl ewood Cliff s , NJ, 1982; 
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J. Pearl. Heuristics : Intelligent Search Strategies for 
Computer Problem Solving. Addison-Wesley , Reading, MA, 1984; 
C.R. Reeves ed. Modern Heuristic Techniques for Combinatorial 
Problems. Halstead Press, Wiley, NY, 1993) . For the three- 
dimensional assignment problem, Pierskalla (W. Pierskalla. The 
tri- substitution method for three-dimensional assignment 
problem. Journal du CORS f 5:71-81, 1967) developed the tri- 
substitution method, which is a variant of the simplex method. 
Frieze and Yadegar (A.M. Frieze and J. Yadegar. An algorithm 
for solving 3 -dimensional assignment problems with application 
to scheduling a teaching practice. Journal of the Operational 
Research Society, 39:989-955, 1981) introduced a method based 
on Lagrangian relaxation in which a feasible solution is 
recovered using information provided by the relaxed solution. 
Our choice of approaches is strongly influenced by the need to 
balance real-time performance and solution quality. 
Lagrangian relaxation based methods have been used 
successfully in prior tracking applications (K.R. Pattipati, 
S. Deb, and Y. Bar-Shalom. A s-dimensional assignment 
algorithm for track initiation. In Proceedings of the IEEE 
Systems Conference, Kobe, Japan, pages 127-130, 1992; Y. Bar- 
Shalom, S. Deb, K.R. Pattipati, and H. Tsanakis. A new 
algorithm for the generalized multidimensional assignment 
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problem. In Proceedings of the IEEE International Conference 
on Systems, Math, and Cybernetics, Chicago, pages 132-13 6, 
1992; A.B. Poore and N. Rijavec. A numerical study of some 
data association problems arising in multitarget tracking. 
5 Large Scale Optimization: State of the Art, W.W. Eager, D.W. 
Hearn and P.M. Pardalos, editors. Kluwer Academic Publishers 
B.V., Boston, pages 339-361, 1994; A.B. Poore and N. Rijavec. 
Partitioning multiple data sets: multidimensional assignments 
and lagrangian relaxation. In P.M. Pardalos and H. Wolkowicz, 

10 editors, Quadratic assignment and related problems: DIMACS 
Series in Discrete Mathematics and Theoretical Computer 
Science, volume 16, pages 25-37, 1994; A.B. Poore and N . 
Rijavec. A lagrangian relaxation algorithm for 
multidimensional assignment problems from multitarget 

15 tracking. SIAM Journal of Optimization, 3 , No. 3:544-563, 
1993) . An advantage of these methods is that they provide 
both an upper and lower bound on the optimal solution, which 
can then be used to measure the solution quality. These works 
extend the method of Frieze and Yadegar (A.M. Frieze and 

2 0 J. Yadegar. An algorithm for solving 3 -dimensional assignment 
problems with application to scheduling a teaching practice. 
Journal of the Operational Research Society, 32:989-995, 1981) 
to the multidimensional case. 
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IV. 2. Probabilistic Framework for Data Association. 



The goal of this section is to explain the formulation of 
the data association problems that governs large classes of 
data association problems in centralized or hybrid 
centralized-sensor level multisensor/multitarget tracking. 
The presentation is brief; technical details are presented for 
both track initiation and maintenance in the work of (A.B. 
Poore. Multidimensional assignments and multitarget tracking, 
partitioning data sets, I.J. Cox, P. Hansen, and B. Julesz, 
editors, DIMACS Series in Discrete Mathematics and Theoretical 
Computer Science, American Mathematical Society, Providence, 
R.I., 19:169-198, 1995; A.B . Poore. Multidimensional 
assignment formulation of data association problems arising 
from multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications, 3:27-57, 1994). 
The formulation is of sufficient generality to cover the MHT 
work of Reid, Blackman and Stein (S.S. Blackman. Multiple 
Target Tracking with Radar Applications. Artech House, 
Norwood, MA, 1986) and modifications by Kurien (T.Kurien. 
Issues in the designing of practical multitarget tracking 
algorithms. In Multi tar get -Multi sensor Tracking; Advanced 
Applications by Y. Bar-Shalom. Artech House, MA, 1990) to 
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include maneuvering targets. As suggested by Blackman (S.S. 
Blackman. Multiple Target Tracking with Radar Applications . 
Artech House, Norwood, MA, 1986) , this formulation can also be 
modified to include target features (e.g., size and type) into 
the scoring function. The recent work (A.B. Poore and O.E. 
Drummond. Track initiation and maintenance using 
multidimensional assignment problems. In D.W. Hearn, W.W. 
Hager, and P.M. Pardalos, editors, Network Optimization, 
volume 450, pages 407-422, Boston, 1996. Kluwer Academic 
Publishers B.V.) significantly extends the work of this 
section to new approaches for multiple sensor centralized 
tracking. Future work will involve extensions to track-to- 
track correlation. 

The data association problems for multisensor and 
multitarget tracking considered in this work are generally 
posed as that of maximizing the posterior probability of the 
surveillance region (given the data) according to 



where Z N represents N data sets, y is a partition of indices 
of the data (and thus induces a partition of the data) , r* is 
the finite collection of all such partitions, r is a discrete 

153 



Maximize 




(4.2) 



CSUR.01USRI 



random element defined on r*, and P(F = y\z N ) is the posterior 
probability of a partition y being true given the data Z N . The 
term partition is defined below; however, this framework is 
currently sufficiently general to cover set packings and 
coverings (A.B. Poore . Multidimensional assignment formulation 
of data association problems arising from multitarget tracking 
and multisensor data fusion. Computational Optimization and 
Applications, 3:27-57, 1994; A.B. Poore. Multidimensional 
assignments and multitarget tracking, partitioning data sets, 
I.J. cox, P. Hansen, and B. Julesz, editors, DIMACS Series in 
Discrete Mathematics and Theoretical Computer Science, 
American Mathematical Society, Providence, R.I., 19:169-198, 
1995) . 

Consider N data sets Z(k) (Jr=l,...,N) each of M k reports 



respectively. In multisensor data fusion and multitarget 
tracking the data sets Z(k) may represent different classes of 
objects, and each data set can arise from different sensors. 
For track initiation the objects are measurements that must be 
partitioned into tracks and false alarms. In our formulation 




and let Z N denote the cumulative data set defined by 




(4.3) 
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of track maintenance (A.B. Poore. Multidimensional assignment 
formulation of data association problems arising from 
multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications, 3:27-57, 1994; 
5 A.B. Poore. Multidimensional assignments and multitarget 
tracking, partitioning data sets, I.J. cox, P. Hansen, and B. 
Julesz, editors, DIMACS Series in Discrete Mathematics and 
Theoretical Computer Science, American Mathematical Society, 
Providence, R.I., 19:169-198, 1995), which uses a moving 

10 window over time, one data set will be tracks and remaining 
data sets will be measurements which are assigned to existing 
tracks, as false measurements, or to initiating tracks. In 
sensor level tracking, the objects to be fused are tracks 
(S.S. Blackman. Multiple Target Tracking with Radar 

15 Applications. Artech House, Norwood, MA, 1986) . In 
centralized fusion (S.S. Blackman. Multiple Target Tracking 
with Radar Applications. Artech House, Norwood, MA, 1986), the 
objects may all be measurements that represent targets or 
false reports, and the problem is to determine which 

2 0 measurements emanate from a common source. 

We specialize the problem to the case of set partitioning 
(A.B. Poore. Multidimensional assignment formulation of data 
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association problems arising from multitarget tracking and 
multisensor data fusion. Computational Optimization and 
Applications, 3:27-57, 1994; A.B. Poore . Multidimensional 
assignments and multitarget tracking, partitioning data sets, 
5 I.J. cox; P. Hansen, and B. Julesz, editors, DIMACS Series in 
Discrete Mathematics and Theoretical Computer Science, 
American Mathematical Society, Providence, R.I., 19:169-198, 
1995) defined in the following way. First, for notational 
convenience in representing tracks, we add a zero index to 
10 each of the index sets in Equation (4.3) , a dummy report z% to 
each of the data sets Z(k) in Equation (4.3), and define a 



value of 0. A partition of the data will refer to a collection 
of tracks of data wherein each report occurs exactly once in 

15 one of the tracks of data and such that all data is used up; 
the occurrence of dummy report is unrestricted. The dummy 
report z 0 serves several purposes in the representation of 
missing data, false reports, initiating tracks, and 
terminating tracks (A.B. Poore. Multidimensional assignment 

20 formulation of data association problems arising from 
multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications, 3:21-51, 1994; 
A.B. Poore. Multidimensional assignments and multitarget 



"track of data 



as 




where 1 k can now assume the 
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tracking, partitioning data sets, I.J. cox, P. Hansen, and B. 
Julesz, editors, DIMACS Series in Discrete Mathematics and 
Theoretical Computer Science, American Mathematical Society, 
Providence, R.I., 19:169-198, 1995). The reference partition 
is that in which all reports are declared to be false. 

Next under appropriate independence assumptions one can 
show (A.B. Poore. Multidimensional assignment formulation of 
data association problems arising from multitarget tracking 
and multisensor data fusion. Computational Optimization and 
Applications, 3:21-57, 1994; A.B. Poore. Multidimensional 
assignments and multitarget tracking, partitioning data sets, 
I.J. cox, P. Hansen, and B. Julesz, editors, DIMACS Series in 
Discrete Mathematics and Theoretical Computer Science, 
American Mathematical Society, Providence, R.I., 19:169-198, 
1995) 

P(r=y|Z N ) T n 

— — — = L = L (4 4) 

P(r- Y °|Z N ) Y (V 2 *> 6 Y VV 

where y° is a reference partition, L i v ..i N is a likelihood ratio 

containing probabilities for detection, maneuvers, and 

termination as well as probability density functions for 

measurement errors, track initiation and termination. Then if 
c. . = -InL. . 
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•In 



P(T=y|Z n ) 

p(r=Y°|z N ) 



(4.5) 



Using Equation (4.4) and the zero-one variable 

Z i l -± N 1 ^i'"*'^ 6 y and 0 otherwise, one can then write 

the problem (4.4) as the following Itf-dimensional assignment 
problem: 

Minimize }] c. . z. . 

Subject To Y, z ii =1 (i^l,...,^), (4.6) 

X z i...i =1 U=l r -,M and p = 2,...,tf-l), 

X l-- L p-1- L p+1 

- L l i 2- :L W-1 

V, 610 ' 11 all i,,...,^, 



where c 0 ... 0 is arbitrarily defined to be zero. Here, each group 
of sums in the constraints represents the fact that each-non- 
dummy report occurs exactly once in a "track of data." One 
can modify this formulation to include multi -assignments of 
one, some, or all the actual reports. The assignment problem 
Equation (4.6) is changed accordingly. For example, if is 
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to be assigned no more than, exactly, or no less than 
n± k times, then the w = l" in the constraint (4.6) is changed to 
< r - r >n ± * respectively. Modifications for group tracking 

and multi -resolution features of the surveillance region will 
5 be addressed in future work. In making these changes, one 
must pay careful attention to the independence assumptions 
that need not be valid in many applications. 

Expressions for the likelihood ratios L ± lf ~.,i N can be found 
in the work of Poore (A.B. Poore. Multidimensional assignment 

10 formulation of data association problems arising from 
multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications, 3:27-57, 1994; 
A.B. Poore. Multidimensional assignments and multitarget 
tracking, partitioning data sets, I.J. cox, P. Hansen, and B. 

15 Julesz, editors, DIMACS Series in Discrete Mathematics and 
Theoretical Computer Science, American Mathematical Society, 
Providence, R.I., 19:169-198, 1995). These expressions 
include the developments of Reid, Kurien (T. Kurien. Issues in 
the designing of practical multitarget tracking algorithms. In 

2 0 Multi target -Multi sensor Tracking: Advanced Applications by Y. 
Bar-Shalom. Artech House, MA, 1990), and Stein and Blackman 
(S . S . Blackman . Mul tiple Target Tracking wi th Radar 
Applications. Artech House, Norwood, MA, 1986). What's more, 
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they are easily modified to include target features and to 
account for different sensor types. In track initiation, the 
N data sets all represent reports from N sensors, possibly all 
the same. For track maintenance , we use a sliding window of 
5 N data sets and one data set containing established tracks. 
(A.B. Poore. Multidimensional assignment formulation of data 
association problems arising from multitarget tracking and 
multisensor data fusion. Computational Optimization and 
Applications, 3:27-57, 1994; A.B. Poore. Multidimensional 

10 assignments and multitarget tracking, partitioning data sets, 
I.J. cox, P. Hansen, and B. Julesz, editors, DIMACS Series in 
Discrete Mathematics and Theoretical Computer Science, 
American Mathematical Society, Providence, R.I., 19:169-198, 
1995) . The formulation is the same as above except that the 

15 dimension of the assignment problem is now N+l . 

IV. 3. Overview of the Lagrangian Relaxation Algorithm. 

Having discussed the N-dimensional assignment problem 
(3.1), we now turn to a description of the Lagrangian 
20 relaxation algorithm. The algorithm will proceed iteratively 
for a loop k=l,. . ., N-2 . At the completion, there remains 
one two-dimensional assignment problem that provides the last 
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step which yields an optimal (sometimes) or near-optimal 
solution to the original iV-dimensional assignment problem. In 
step k of this loop (summarized in Section IV. 3) one starts 
with the following (N-k+1) -dimensional assignment problem with 
one change in notation. If then the index notation 1 2 and 

L 1 are to be replaced by i 1 and M lf respectively. 



Minimize £ cJT 1 ± z 2 H ? +1 ; 

7 . . • L k' L k+l ' * * X W ■ L k J -k+l ' * • J ~N 

1 k 1 k+l' * ' 1 N 

Subject To E . <C---i N =1 ' i * = 1 "-- I V ( 4 - 7 ) 

E N ~ k+1 -I ' -1 M 

. Zl k i k+l i k+2'- ± N~ ' 1 k + l~ 1 ' " • " ' W Jt + l' 
± k X k^ ' * * 1 W 



for i = 1, . . . ,M and p = £+2, . . . ,iV-l, 



To ensure that a feasible solution of Equation (4.7) always 
exists for a sparse problem, all variables with exactly one 
nonzero index (i.e., variables of the form zj^ 1 for 
l k =l, . . . ,L k and z 0 M ..* + o.„ 0 for i p =l, . . . , m . and p=Jc+l, . . . ,N) are 
assumed free to be assigned with the corresponding cost 
coefficients being well-defined. This assumption is valid in 
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the tracking environment (A.B. Poore. Multidimensional 
assignments and multitarget tracking, partitioning data sets, 
I.J. cox, P. Hansen, and B. Julesz, editors, DIMACS Series in 
Discrete Mathematics and Theoretical Computer Science f 
5 American Mathematical Society, Providence, R.I., 19:169-198, 
1995; A.B. Poore. Multidimensional assignment formulation of 
data association problems arising from multitarget tracking 
and raulti sensor data fusion. Computational Optimization and 
Applications, 3:27-57, 1994). 

10 Section IV. 3.1 presents some of the properties associated 

with the Lagrangian relaxation of (4.7) based on relaxing the 
last (N-Jc) -sets of constraints to a two-dimensional one. 
Section IV. 3. 2 describes a new approach to the problem of 
recovering a high quality feasible solution of the original 

15 (N-k+1) -dimensional problem given a feasible solution (optimal 
or suboptimal) of the relaxed two-dimensional problem is 
described hereinafter. A summary of the relaxation algorithm 
is given in Section IV. 3. 3, and in Section IV. 3. 4, we 
establish the maximization of the Lagrangian dual (an 

2 0 important aspect of the relaxation procedure) to be an 
unconstrained nonsmooth optimization problem and then present 
a method for computing the subgradients . 
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IV. 3.1 The Lagrancrian Relaxed Assignment Problem 
The (N-k+1) -dimensional problem (4.7) has {N~k+1) sets 
of constraints. A (Af p +1) -dimensional multiplier vector (i.e., 
u p G R m p +1 ) associated with the p th constraint set will be 
denoted by u p = (ug, uf , . . . , u%) T with u% = 0 being fixed for each 
p = k + 2, . . . ,N and included for notational convenience only. 
Now, the -dimensional assignment problem (4.7) is 

relaxed to a two-dimensional assignment problem by 
incorporating {N-k-1) sets of constraints into the objective 
function via the Lagrangian. Although any {N-k-1) constraint 
sets can be relaxed, we choose the last {N-k-1) sets of 
constraints for convenience. The relaxed problem is 
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& N - k+1 (u k+2 , ■ ■ - r u N ) s Minimize <P N _ k+1 ( z w "* +1 ; u k+2 , . . . u 
Minimize £ cf** 1 ., z.T-i 



p=£+2 i p =0 



,P 



Minimize 



p=*+2 i p =0 



N-k+1 - 



- L Jt x Jt+i ^-w 



= 0 

p=*-2 z p =0 

Subject To X) zr/ +1 i i = If 1,= 1,...L„ 



1 k+l :L k+2'" X N 



7 . - J -Jc- L Jt+l i Jt+2 N ' £+1 ' ' /C+l 



(4,8 



One of the major steps in the algorithm is the maximization of 
0 N-k+i ( uk+2 r • • u N ) with respect to the multipliers 
[u k+2 f . . . , u N ) . It turns out that & N „ k+1 is a concave, 
continuous, and piecewise affine function of the 
multipliers (u k+2 r . . u N ) , so that the maximization of & T , , is 
a problem of nonsmooth optimization. Since many of these 
algorithms require a function value and a subgradient of 0 , , 

3 N-k+l 

at any required multiplier value (u* +2 , . . . , u w ) , we address 
this problem in the next subsection. We note, however, that 
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there are other ways to maximize @ N _ k+1 and the next subsection 
just addresses one such method. 



IV. 3. 2 Properties of the Lagrangian Relaxed Assignment Problem 
5 For a function evaluation of &„ . , , we show that an 

N-k+l ' 

optimal (or suboptimal) solution of this relaxed problem (4.8) 
can be constructed from that of a two-dimensional assignment 
problem. Then, the nonsmooth characteristics of $ N „ k+1 are 
addressed, followed by a method for computing the function 
10 value and a subgradient . 

Evaluation of @ N . k+1 * Define for each (l k/ i k+1 )an index 
j N ) =Uk +2 (lkf ik + i) f»,j N (lk'ik + i) ) and a new cost function 



(j*+2 ^k' 1 k + l) ' * * * ' ^N^k' ^Jfc+l^ ^ " 

hT\ i + £ ufn = o,i,...,m ; 

[ and p=k+2, . . . ,N 



argrmim ° 7 



14.9) 



2 W-Jfc + 1 

^k^k+1 ^k^k^k+2 ^k' * * ' Jff^Jf ^Jt+l^ 



+ £ ^d„i M ) for (1,,!^)* (0,0), 
= X) Minimum <0,c i + S u f k 
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Given an index pair {l k , i k+1 ) , (j k+2 , . . . , j N ) need not be unique, 
resulting in the potential generation of several subgradients 
(4 . 17) . Then, 

S (u k * 2 , ,u N ) = 

N-k+l K ' ' " * ' U ' 

Minimize <P N _ k+1 (z 2 ; u k+2 , . . . , u N ) 



L k U 

z l^ e{0 ' 1} for a11 1 k' i k+l- 



Subject To 22 z i± ==1 f 1 iT 1 f 



-*■/:+ 1 U 



(4.io; 



As an aside, two observations are in order. The first is that 
5 the search procedure needed for the computation of the relaxed 
cost coefficients in (4.9) is the most computationally 
intensive part of the entire relaxation algorithm. The second 
is that a feasible solution z N ~ k+1 of a sparse problem (4.7) 
yields a feasible solution z 2 of (4.10) via the construction 



10 



2 



1, it z 1 ± i m± = 1 for some »■ , j.„) , 

0, otherwise. 



thus, there are generally solutions other than the one nonzero 
index solution. 
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The following Theorem 4 .1 gives a method for evaluating @ N _ k+1 
and states that an optimal solution of (4.8) can be computed 
from that of (4.10) . Furthermore, if the solution of either of 
these two problems is E-optimal, then so is the other. The 
5 converse is contained in Theorem 4.2. 

Theorem 4.1. Let w 2 be a feasible solution to the two- 
dimensional assignment problem (4.10) and define w^"^ 1 by 

N-k+l 2 ' -f { ' - \ _ / ■ • \ 

and (l k ,i k+1 )*(0,0; 

™TJili k+z ..N= 0 if (i k+2 , . . . i N ) * (j k+2 , j N ) 

and (l k ,i k+1 )*(0,0: 

N-k+l - r N-k+l , \^ p ~ 



N-k+l A . ^ N-/C+1 , \^ P ^ r\ 

W 00i «.i = 0 lf C 00i .-I + Zy Ui > 0. 



(4.11) 



Then w N ~* :+I is a feasible solution of the Lagrangian relaxed 
problem and 

1 0 <P N . k+1 ( w*-* +1 / u k+2 , ... , u N ) = <p- k -! ( n 2 ; u k+2 , . . . , u N ) . 

If, in addition, w 2 is optimal for the two-dimensional problem, 
then w*^* 1 is an optimal solution of the relaxed problem and 
®N- k+ Au k + 2 , ...,u N )= S_ k+1 (u k + 2 , ...,u N ). 
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Proof. Let w*"** 1 and w 2 be as in the hypotheses , and let 

<P*- k +iW' k+1 ;u** 2 , • • wu") and ' 

denote the objective function values of (4.8) and (4.10) 7 
respectively. Direct verification shows that w®~ k+1 satisfies 
the constraints in (4.8) and 

<P*- k+ i (w^/u** 2 ; ...,1/0= (w 2 /^ 2 , ...,11") . 

For the remainder of the proof, assume that w 2 is optimal for 
(4.10). Let x^"^ 1 satisfy the constraints in (4.8) and define 

x ii = Xi i x iT\ -i for (0,0) , 

x 0 2 0 = l if (l k ,i k+1 ) = (0,0) and 

£, K s0 forsome (i — v. (4 - 12) 

x 0 2 0 = 0 if (J^i ) = (0,0) and 



^'-■i, + £ u£ >0 for all U, +2 ,-, i M ) . 



Note that x 2 satisfies the constraints in (4.10) and 
fcr-jc+i (x™;^ 2 , . . . , u w ) > (x 2 ;^ 2 , . . . , u N ) 

= <P N - k+ A^- k+1 ;u k+2 f ...,u N ) 
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for any feasible solution xp~ k+1 of the constraints (4.8) . This 
implies v^~ kJrl is an optimal solution of (4.8). Next, 

(u* +2 , - • • / vP) = S N _ k+1 (u k+2 , ...,u N ) 
follows immediately from 
5 0 N - k+1 (w*-** 1 ;^ 2 , . . . , u N ) = <f N _ k+1 {^;u k+2 , . . . , vP) 

for an optimal solution w 2 of (4.10) . With the exception of one 
equality being converted to an inequality, the following 
theorem is a converse of Theorem 4.1. 

Theorem 4.2. Let w^~ k+1 be a feasible solution to problem (4.8) 
10 and define w 2 by 

™li = £ w iT\ ...i for (l.,i. +1 )*(0,0) , 

■ L k x k*l . . 1 k 1 k*i 1 k+2 ^-W k' k+1' ' ' 

1 k*2 "' 1 U 

Wq 0 = 1 if [l k ,i k+1 ) = (0,0) and 

c ooi M .-i„ + 2, ^0 some (i^,-, ± N ) , 

15 p= * +2 

^ 0 2 0 = 0 if ll k ,i k+1 ) = (0,0) and 

Then w 2 a 
feasible 

2 0 solution of the problem (4.10) and 

0N- k+ l (^ k+1 ;u k + 2 , ...,12*) > <?Wi (w 2 /^ 2 , . . . , u N ) . 
If, in addition, w^- k+1 is optimal for (4.8), then w 2 is an 
optimal solution of (4.10), 
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^- Jt+ i(w ff -* +1 /u k+2 / ... / u tf ) = 4-* + i(^;u* +2 , ...,u*) and 

(^ +2 / . . . / u N ) = S N _ k+1 (u k+2 , ...,u N ). 
Proof: Let vJ*~ k+1 and w 2 be as in the hypotheses, and let 

fa-jc+iiv*'** 1 /^ 2 , • • • / u") and (wW 2 , ...,11") denote the 

5 objective function values of (4.8) and (4.10), respectively. 
Direct verification shows that w 2 satisfies the constraints in 
(4.10) and 

0*-*+i (w^/u** 2 , . . . , u") > 0 N _ k+1 ( w 2 /^* 2 , . . . , u N ) . 
For the remainder of the proof, assume that w* - ** 1 is optimal 
10 for (4.8) and construct w 2 as above. Using Theorem 4.1, 
construct w N ~ k+1 by replacing w*' k+1 in (4.11) with w N ~ k+1 . Then, 
from that theorem 

0 N - k+ i ( * N ~ k+1 ;u k + 2 , . . . , u N ) = ^ k+1 {w 2 ;^ 2 , ...,u N ) 

15 Optimality of vF'** 1 then implies the last inequality is in fact 

an equality. This proves the theorem. 

The Nonsmooth Optimization Problem. Next, we address the 

nonsmooth properties of the function <P N „ k+1 ( u k+2 , . . . , u*) as 

explained in the following theorem. 
20 Theorem 4.3. (G.L. Nemhauser and L.A. Wolsey. Integer and 

Combinatorial Optimization, Section II. 3. 6. Wiley, New York, 

NY, 1988). Let <P N _ k+1 be as defined in (4.7), let V N _ k+1 (z N ' k+1 ) 

170 



CSUR.01USRI 

be the objective function value of the (N-k+1) -dimensional 
assignment problem in equation (4.7) corresponding to a 
feasible solution z N ' k+1 of the constraints in (4.7), and let 
£N-k + i be an opt i ma i solution of (4.7). Then, <P N _ k+1 (u k+2 , . . . , u N ) 
5 is piecewise affine, concave and continuous in {u k+2 , . . .,u N J and 

t> N _ k+1 (u k +\ . . . , u w ) * V N _ k+1 (z^) < V N _ k+1 (z^) . (4 . 14 ) 

Most of the algorithms for nonsmooth optimization are based on 
10 generalized gradients, called subgradients , given by the 
following definition for a concave function. 

Definition 4.1. At D = (u k+2 , . . . , u N ) the set d$ N _ k+1 {u) is called 
a subdifferential of @ N _ k+1 and is defined for the concave 
15 function & N _ k+1 (u) by 

d# N _ k+1 (u) ={ge R M ^ +1 x ... x R M « +1 | {w) ~0 N _ k+1 (u)< (4.15) 

g T {w-u) for all w e R Mk+2+l x.»xR^ +1 }, 

where = = = 0 are a11 permanently fixed. (Recall that 
20 these were used for notational conveience only.) A vector 
g e d@ N _ k+1 is called a subgradient. 
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There is a large literature on such problems, e.g., 
(J.-B. Hiriart-Urruty and C. Lemarechal . Concex Analysis and 
Minimization Algorithms I& II. Springer-Verlag, Berlin, 1993; 
K.C. Kiwiel. Methods of descent for nondif f erentiable 
5 optimization. In Lecture Notes in Mathematics 1133, A. Bold 
and B. Eckmann, eds. Springer-Verlag, Berlin, 1985; 
C. Lemarechal and R. Mifflin. Nonsmooth Optimization. Pergamon 
Press, Oxford, UK, 1978; B.T. Polyak. Subgradient method: 
A survey of Soviet research. In C. Lemarechal and R. Mifflin, 

10 editors, Nonsmooth Optimization, pages 5-29, N.Y., 1978. 
Pergamon Press.; H.Schramm and J. Zowe. A version of the 
bundle idea for minimizing a nonsmooth function: Conceptual 
idea, convergence analysis, numerical results. SI AM Journal on 
Optimization, 2, No. 1 : 121-152 , February, 1992; N.Z. Shor. 

15 Minimization Methods for Non-Diff erentiable Functions. 
Springer-Verlag, Nem York, 1985; P.Wolfe. A method of 
conjugate subgradients for minimizing nondif f erentiable 
functions. Mathematical Programming Study, 3:145-173, 1975), 
and we have tried a variety of methods including subgradient 

2 0 methods (N.Z. Shor. Minimization Methods for Non- 
Diff erentiable Functions. Springer-Verlag, New York, 1985) 
and bundle methods (J.-B. Hiriart-Urruty and C. Lemarechal. 
Convex Analysis and Minimization Algorithms I& II. Springer- 
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Verlag, Berlin, 1993; K.C. Kiwiel . Methods of descent for 
nondif ferentiable optimization. In Lecture Notes in 
Mathematics 1133, A. Dold and B.Eckmann, eds . Springer -Verlag, 
Berlin, 1985) . Of these, we have determined that for a fixed 
number of nonsmooth iterations (e.g., twenty), the bundle 
trust method of Schramm and Zowe (H.Schramm and J. Zowe. A 
version of the bundle idea for minimizing a nonsmooth 
function: Conceptual idea, convergence analysis, numerical 
results. SIAM Journal on Optimization, 2, No. 1:121-152, 
February, 1992) provides excellent quality solutions with the 
fewest number of function and subgradient evaluations, and is 
therefore our currently recommended method. 

An Algorithm for Computing & N _ k+1 a.nd a Subgradient. 
Most current software for maximizing the concave function 
1 re quires the value of the function and a subgradient at 
a point {u k+2 , . . . , u N ) . Based on the previous two sections, 
this can be summarized as follows. 

1. Starting with problem (4.7), form the relaxed 
problem (4.8) ; 

2. To solve (4.8) optimally, defined the two- 
dimensional assignment problem (4.10) via the transformation 
(4.9) ; 

3. Solve the two-dimensional problem (4.10) optimally; 
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4. Reconstruct the optimal solution, say w 
(4.8) via equation (4.9) as in Theorem 4.1; 

5. Compute the function 



EN-k+1 *N-k+l 



1 k j -k+l r '"j-N 



p=/c+2 i p =0 



w l i -i " 



(4.16) 



10 



6. A subgradient is given by 



9£ 



/c+2 



• i jt 1 Jt+l" - ^p-l^p+l"" *n 



i -i "I 



for i p = 1/ . . . , M and p = £+2, 



(4-17) 



Several remarks are in order. First, the optimal solution of 
the minimization problem (4.8) is required before one can 
remove the minimum sign, replace z N ~ k+1 by the minimizer $®~ k+1 
and differentiate with respect to u[ to obtain a subgradient 
15 as in (4.17). If vF~ k+1 is an approximate solution of (4.8), 
then the subgradient and function values are only approximate 
with accuracy depending on that of ^' k+1 . Although one can 
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evaluate the sums in (4.16) and (4.17) in a straight forward 
manner, another method is based on the following observation. 
Given the multiplier vector ( u k+2 r * * * ' ' 
Ul k {l k+l ) , i k+1 {l k+1 ) )Y k+l be an enumeration of indices |l k /i k+1 } of 
5 w 2 (or the first 2 indices of w N ~ k+1 constructed in equation 
(4.9)) such that < |lwUwllt/ l for (l k d k+1 ) , d k+1 ) ) * (0, 0) 
and (l k {l k+1 ) , i k+1 (l k+1 ) ) = (0,0) for l t+1 =0 regardless of whether 

2 

p\r 00 -l or not. (The latter can only improve the recovered 
feasible solution.) The evaluation of the bracketed quantity 
10 in (4.17) for a specific i p ^l and p = k+2 , . . . ,N is one minus 
the number of times the index value i p appears in the (p-k+l) th 
position of the (N-k+1) -tuple in the list 

k^k+l^ ' ^Jfc+l ^Jt+l^ r 3k+2' ' ' ' '^w} 7 * +1 _ n • 

Finally, we have presented a method for computing one 

15 subgradient. If v^'^ 1 is the unique optimal solution of (4.8), 
then ® N - k+1 {u) is dif f erentiable at u, g is just the gradient at 
u, and the subdif f erential d& N _ k+1 ( u) = { g) is a singleton. If, 
corresponding to u, $*~ k+1 is not unique, then there are 
finitely many such solutions, say { i^~ k+1 (1) , . . . , v^~ k+1 ( m ) } with 

20 respective subgradients {g(D , - - - ,g(m) } r the subdif f erential 
d&(u) is the convex hull of {g (1) , . . . , g (m) } (J.L. Goffin. On 
convergence rates of subgradient optimization methods, 
mathematical programming. Mathematical Programming, 13:329- 
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347, 1977). These nonunique solutions arise in two ways. 
First, if the optimal solution of the two-dimensional 
assignment problem is not unique, then one can generate all 
optimal solutions of the two-dimensional assignment problem 
5 (4.10) . Corresponding to each of these solutions and to each 
index pair (1 , i k+1 ) in each solution, compute the indices 
( 3 k+2 r - • • t J N ) in (4.9). If these j's are not unique, then we 
can enumerate all the possible optimal solutions ^~ k+1 of 
(4.8). Given these solutions, one can then compute the 
10 corresponding subgradients from (4.17). In most nonsmooth 
optimization algorithms, one uses an e-subdif f erential . 

Definition 4.2. At u = (u k+2 , . . . , u N ) the set d € @ N - k+1 (u) is 
called an e-subdifferential of 0 AT . t and is defined for the 

N-k+l 

15 concave function @ N _ k+1 (u) by 

5 e %- k+ i ( u) = { ff e K M * +2+1 * • • • * 1 <V**i ( " <Wi (u)< 
g T (w-u) + e for all w eM Mfc - 2+1 * . . . x R^ 1 }, 

where = Wq = Uq = 0 are all permanently fixed. (Recall that 
these were used for notational convenicence only.) A vector 
g e d e 0 N _ k+1 (u) is called an e-subgradient. 
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IV. 3. 3 The Recovery Procedure. 
The next objective is to explain a recovery procedure, i.e., 
given a feasible (optimal or suboptimal) solution w 2 of (4.10) 
(or v/*~ k+1 of (4.8) constructed via Theorem 4.1), generate a 
5 feasible solution of equation (4.7) which is close to w 2 

in a sense to be specified. We first assume that no variables 
in (4.7) are preas signed to zero; this assumption will removed 
shortly. The difficulty with the solution w* 7- ** 1 in Theorem 4.1 
is that it need not satisfy the last (N-k-1) sets of 

10 constraints in (4.7) . (Note, however, that if w 2 is an optimal 
solution (4.10) and \J*~ k+1 , as constructed in Theorem 4.1, 
satisfies the relaxed constraints, then w®~ k+1 is optimal for 
(4.7) .) The recovery procedure described here is designed to 
preserve the 0-1 character of the solution w 2 of (4.10) as far 

15 as possible: If w 7 , =1 and l.*0 or i. -*0, the 

corresponding feasible solution z N ~ k+1 of (4.7) is construed so 
that z^ i ^ 1 „ i ^ = l for some (i k+2 , . . . , i N ) . By this reasoning, 
variables of the form z n w n ~* +1 , can be assigned a value of one 
in the recovery problem only if Wq 0 ~1. However, variables 

20 z ooi**«.i will be treated differently in the recovery procedure 
in that they can be assigned either zero or one independent of 
the value w$ 0 . This increases the feasible set of the 
recovery problem, leading to a potentially better solution. 
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Let I {l k (l k+1 ) , ± k+1 (l k+l ) \ Lk+1 be an enumeration of indices 

{l k , i k+1 } of w 2 (or the first 2 indices of w^ _jc+1 constructed in 
equation (4.11)) such that 

<<i, +l) ,i w u, +1 ) = l for 

^ 1 k( 1 k + i^i k+1 il k+l )) = (0 r 0) for 
regardless of whether ^ 0 2 0 = 1 or not. To define the (N-k) - 
dimensional assignment problem that resores feasibility, let' 



^1 i i ~ ^1 f 1 ) i ( 1 ) i i fOJO k — 1 

■ L 2- L 3*'*- L N i ^ 2 ' 2 * 2 ' 3 * * * W 

N-k N-k+1 . A - 0 . 

°i i i = c i a )i a )i i 4.18 

— c • • * * ■ * 

1 1 ^2...Jc+l) X 2 ^2...k+l^ Z 3 ^3...Ar+l^ " ' m 1 k^kk+^) Z k+1 ^Jt+l' 1 k+2 ' ' ' Z N 



for k^2 and I = 0, . . . , L 



where 

^ ...m = V 1 -! 0 • • • 0 ^ ( ) ■ K < (••• < ^ ( ))■••)) (4.i9: 



10 for m - 2 , . . .Jc+l and where denotes function composition. 

For example, 1 23 = 1 2 {1 3 ) and 1 234 = 1 2 °1 3 (1 4 )= 1 2 (1 3 (1 4 )). 

Then the (N- k ) -dimensional assignment problem that 
restores feasibility is 
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Minimize 



E N-k N-k 

, . . J -Jk+l- L Jfc+2* * * X N - L Jt+l- L Jt+2* * '^N 

± k+l J -k+2 * * ' 2 W 

Subject To £ <X 2 ...i N =l. Vr^-'-'V 

EN-k _ -l ' _ -l *, 



(4.20) 



1 k 1 k^l * * ■ 2 p-lVhl * * * 2 W 



W-Jt 



for i p = l/ . - . ,M and p = k+3, . . . ,N-1, 

/ j Z 1 7 . — \ 1 1T = lr , , , M. T | 



W-Jt 



2 i i e {0, 1} for all I. i, . . . , i„. 

Let Y be an optimal or feasible solution to this [N-k) - 
dimensional assignment problem. The recovered feasible 
solution is defined by 



1, if i 1 = i 1 d 2- .. Jt+1 ) / W^.. -Jt+1 ), 

i 3 = i 3( I 2...j t+ i)"-w i k = i k^ 1 kk + i^ ( 4 - 21 ) 

0, otherwise. 



Zi 4 



1^3... l ff 



Said in a different way, the recovered feasible solution z N 
5 corresponding to the multiplier set { u k+2 , . . . , u N } is then 
defined by 
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J i 1 U 2 .*+l* ^2 ^2 ^3 ^3- '"^Jt^JtJt+l^ *k*l ^Jk+l' ^k+2'"^N 



1, if y 2 . .... =1, 
fc 0, otherwise, 



where l m ...k + i is defined in (4.19) and "°" denotes function 
composition. 

The next objective is to remove the assumption that all 
5 cost coefficients are defined and all zero-one variables are 
free to be assigned. We first note that the above 
construction of a reduced dimension assignment problem (4.20) 
is valid as long as all cost coefficients cF~ k are defined and 
all zero-one variables in z N ' k are free to be assigned. 
10 Modifications are necessary for sparse problems. If the cost 
coefficient cf/T 1 ., , 7 w , is well-defined and the zero-one 

± k K ' L k+l 1 •'■k+l K± k+1 } ± k+2" J -N 

variable zf/T 1 *, n w 7 . is free to be assigned to zero or 

one, then define cf~ k i ± = cf'^ 1 , . n w , as in (4.18) with 

zf k i ; being free to be assigned. Otherwise, zf~ k , , is 

15 preassigned to zero or the corresponding arc is not allowed in 

the feasible set of arcs. To ensure that a feasible solution 

exists, we now only need ensure that the variables zf~ k n n are 

- L jt+i u "* u 

free to be assigned for l k+1 =0, l,-,L k+1 with the corresponding 
cost coefficients being well-defined. (Recall that all 
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variables of the form zfo-o 1 anc * z q-q1\-q are free (to be 

k p 

assigned) and all coefficients of the form c^qSq an d 
c o M "oi +1 o..-o are well-defined for l k =0,...,L k and i p =0,...,M p for 
p=k+l,...,N.) This is accomplished as follows: If the cost 
5 coefficient c™ i)ijt+i( ^ +i) o...o is well-defined and ^u^ l )i Jfe+1 {i Jt+1 )o...o 
is free, then define c£* 0 „. 0 = ^t^i^ti^jo.-.o with z £*o.-o bein 9 
free. Otherwise, since all variables of the form z^Q^and 
z oi^]o-o are f ree to be assigned with the corresponding costs 
being well-defined, set c 2 *~* 0 ... 0 = Ci* u^joo-.o + c o^u Jt+1 )o-»o * where the 
10 first term is omitted if Ijc f-ljt+iJ =° an <d the second, if 
ijt+i^Jt+i^ 0 - For 1 k( 1 k + i)=° and ijfc + 2^jk+2^=0/ define Cq~£ = c 0 w ; o * +1 . 



IV. 3. 4. The last step k = N-2 . 
Minimize 22 12 °i i z i ±=VAz 2 ) 

Vi = 0 V s o 

Subject To }2 z\ ' ± = 1, 1^ = 1,...,!^, (4.22) 

6 {0,l} for all 



^2 



The description of the algorithm ends with k = JV+2. The 
15 resulting recovery problem defined in (4.10) with k = N-2 is 
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Let Y be an optimal or feasible solution to this two- 
dimensional assignment problem. The recovered feasible 
solution z N is defined by 



V2 = V 2 (Wi)'VrVi(Vi) and y WN =1 
or if (Vi'V = <°'°) ' 
0, otherwise. 



(4.23) 



Said in a different way, the recovered feasible solution ^ is 
then defined by 



2. 



1, if Y 1 = 1, (4.24) 

■ L N-1 ± N 

0, otherwise. 



where l m ... N -i is defined in (4.19) and w o" denotes function 
10 composition. To complete the description, let 

Ul N ^{l N ) r i N {l N )\ LlJ be an enumeration of indices {l N _ lf i N ) of Y 

i N o 

such that ^. i(V , ijr(lN) =l for (l w . x (1 N ) , i w (l N ) ) * (0, 0) and 
(-Vidtf) ,i w d w ))= (0,0) forl w =0 regardless of whether Y 00 =l or 
not. Then the recovered feasible solution can be written as 



15 



1, if i^i^l^J f i 2 = i 2 (l 2 J ,i 3 = i 3 (l^ N ) 

1 N-l ~ 1 N-l (1 N-1N^ ' i N ~ * N (1 ' 

0, otherwise. 

IV. 3. 5. The Upper and Lower Bounds 
The upper bound on the feasible solution is given by 



(4 .25; 



20 
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v**> = E ^-i^ (4 ' 26) 

and the lower by <P N (u 3 , ...,u N ) for any multiplier value (u 3 , u^) . 
In particular, we have 

5 

4 N iu 3 ,...,u N ) < <P N (u 3 ,...,u N ) £ V N {J N )< V N (2 N ), (4.27) 

where {u 3 ,...,u N ) is any multiplier value, (u 3 ,.../^) is a 
maximizer of & N {u 3 , u N ) , ^ is an optimal solution of the 
10 multidimensional assignment problem (4.7) and z N is the 
recovered feasible solution defined by (4.25). Finally, we 
conclude with the observation that V N {z) =V 2 {Y) where Y is the 
optimal solution of (4.23) used in the construction of z in 
equations (4 . 24) - (4 .26) . 

15 

IV. 3. 6. Reuse of Multipliers 
Since the most computationally expensive part of the 
algorithm occurs in the maximization of @ N -k +1 , the development 
of algorithms that can make use of hot starts or multipliers 
2 0 close to the optimal are fundamentally important for real-time 
speed. The purpose of this section is to demonstrate that the 
multiplier set obtained at stage k>l provide good starting 
values for those obtained at step k+1 . 
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Theorem 4.4. Let (u 3 ,...,u N ) denote a multiplier set obtained in 
the maximization of <P N {u 3 , u N ) . Then this multiplier set 
satisfies the string of inequalities 

0 N (u\ . . .,u N ) <0 N _ x (u\ . . .,u N ) <: . . . <0 3 {u N ) <0 2 <V N {z) , (4.28) 

where after the first step in the maximization of <P N the 
multipliers are not changed in the remaining steps. 
Furthermore, to improve this inequiality, let 
[u 2 + k l ^k^i f ^ u N l N^i ) denote a maximizer of $ N _ k+1 (u k+2 , u N ) . 

Then, we have 

0 N _ k+1 (u k+2 , . .. r u N )^ N _ k+1 {u 2+k ' N - k+1 f . .. r u N - N ' k+1 ) (4.29) 

z$ N _ k (u 3+k ' N - k+1 ,...,u N > N - k * 1 ) 
z<P N _ k (u 3 + k ' N - k , . ..,u N ' N - k ) 
<: . . . <0 3 {u N - 3 ) <0 2 <V N (z) . 



Proof: Suppose we have a value of the multipliers {u k+2 , u N ) 
obtained in the maximization of 

$ N - k+1 (u k+2 , ...,u N ) = Minimize 0 N _ k+1 {z N ~ k+1 } u k+2 , . . . , u N ) (4.30) 

Subject To £ z iT\ i =1 > 1 k = 1 ' • • - L k' 

1 k*\ 1 k*f • • 1 H 



EB-k+l =1 i =1 M 



where 
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EN-k+1 



N-k+1 



N ^k^k+l * 



N-k+1 - 
Z l i i ~ 1 



(4.31) 



w-Jt+i ^ p w-Jt+i ^ V s p 



p=k+2 



p=k+2 i p =0 



These need not be the maximizer; however, we do assume that we 
have solved the above minimization problem optimally to 
evaluate @ N _ k+1 ( u k+2 , u N ) . Just as in the definition of the 
recovery problem discussed earlier, let 

5 {^ic^it+i)' ^Jt+i ^Jt+i^ ^1**1=0 be an enumeration of indices {l k ,i k+1 } 
of the optimal solution w 2 of problem (4.10) (or the first 2 
indices of the solution w^~ k+1 of the relaxed problem (4.8) 
constructed in equation (4.11)) such that wf n * n , =1 for 
^(^i^^l^iiH 0 ' 0 ) and (l t (l w ),i M (2 w ))=(0 ( 0) for l k+1 =0 
10 regardless of whether Wq Q = 1 or not. 

Then, the final value of the objective function can be 
expressed as 
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Minimize <p N _ k+1 ( z N ' k+1 ; u k+2 , u N ) (4.32) 
Subject to E. ^U^d.^W-i/ 1 ' = 1 > -> Vi' 



where 



0 n . k+1 (z N - k+1 ;u k+2 ,...,U N ) = (4.33) 



£ £ -£ 



p=*+2 l p =0 



If now another constrain set, say the {k+2) th set, is added to 
this problem, one has 

^jt + i< u * +2 '- uW > 5 Minimize f N _ k+1 (z N ~ k+1 ; u k+2 , u w ) = (4.34) 
Minimize £ _ | cf u J >i (1 + J uf kjj^i^!^)..^ 

- i t < 

p=Jc+2 i p =0 p 

Subject to £ 'iXiMiWWi/ 1 ' ^r 1 '-^' 

2 Jt+2"" 2 W 

Ez w ~* +1 =1 i =1 M 
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Since the constraint /J zf K 1 )± {1 )± ... ± =1 for 

1 k+l 1 k+3" 1 N 

i ^+2 ==1 ' -' M k+2 is now imposed, the feasible region is smaller, so 
that one has 

where the fact that @ N ^ k+1 does not depend on the multiplier set 
u k+2 due the added constraint set. Also, the last equivalence 
follows from the fact that @ Nam jc +1 {u k+3 , ... , u N ) is the relaxed 
10 problem (with k replaced by Jc+1) for the recovery problem in 
step k of the above algorithm. Continuing this way, one can 
compute {u 3 ,...,u N ) at the first step (Jc=l) , fix them thereafter, 
and perform no optimization at the subsequent steps to obtain 



a> M (u 3 / ... / u ir )^ M (uS... / u N ) z~.<<P 3 (u N ) <<P 2 =V N (2) , 
where <P 2 i s defined in (4.22) . 

To explain how to improve this inequality, let 
(u 2+ *' w ~* + V..,u w ' w -* +1 ) denote a maximizer of & N _ k+1 ( u k+2 , u N ) . 
2 0 Then by the same reasoning one has the result as stated in the 
theorem. 
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V. SUMMARY OF THE LAGRANGIAN RELAXATION PROBLEM 



Given the multidimensional assignment problem 



Minimize 




Subject To 




(5.1) 



2 1* * ^p-l^p+l* ' %1 N 




M p and p = 2 



r • • • r 



N-l) r 



Ez. . = 1, (i = 1, . . .MJ , 

7 7 IN 



z. .6(0,1} for all i-, . . . i 



where c 0 „ 0 is arbitrarily defined to be zero and is included 
for notational convenience and where the superscript N on both 
5c and z is not an exponent, but denotes the dimension of the 
subscripts and the assignment problem as stated in (3.5) . In 
relaxed and recovery problems c<£. 0 need not be zero! In this 
problem, we assume that all zero-one variables zf ' m± with 



precisely one nonzero index are free to be assigned and that 
10 the corresponding cost coefficients are well-defined. (This 
is a valid assumption in the tracking environment (A.B. Poore. 
Multidimensional assignments and multitarget tracking, 
partitioning data sets, I.J. Cox, P. Hansen, and B. Julesz, 
editors, DIMACS Series in Discrete Mathematics and Theoretical 
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Computer Science, American Mathematical Society, Providence, 
R.I., 19:169-198, 1995; A.B. Poore . Multidimensional 
assignment formulation of data association problems arising 
from multitarget tracking and multisensor data fusion. 
5 Computational Optimization and Applications, 3:27-57, 1994).) 
Although not required, these cost coefficients with exactly 
one nonzero index can be translated to zero by cost shifting 
(A.B. Poore and N. Rijavec. A lagrangian relaxation algorithm 
for multidimensional assignment problems arising from 

10 multitarget tracking. SIAM Journal of Optimization , 3, No. 
3:544-563, 1993) without changing the optimal assignment. 

Having explained many of the relaxation features, it is 
now appropriate to present the Lagrangian relaxation 
algorithm, which iteratively constructs a feasible solution to 

15 the AT- dimensional assignment problem (5.1). 

Algorithm 5.1 (Lagrangian Relaxation Algorithm) To 
construct a high quality feasible solution, denoted by z N , of 
the assignment problem (5.1), proceed as follows: 

0. Initialize the multipliers {u k+2 ,...,u N ) , e.g., 
20 (u k+2 ,... f u N ) = (0,...,0) . 

For k-l,...,N-2, do 

1. For the Lagrangian relaxed problem (4.8) from the problem 
(4.7) by relaxing the last (N-k-1) sets of constraints. 
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2. Use a nonsmooth optimization technique to solve 

Maximize { 0 N „ k+1 1 u p ER Mp+1 for p = k+2,...,N (5.2) 

with u 0 p =0 being fixed) , 

where %_ k+1 ( u k+2 r u N ) is defined by equation (4.8). The 
algorithm in Section IV. 3. 2 provides one method for 
computing a function value and a subgradient out of the 
subdifferential at [u k+2 , u N ) , as required in most 
nonsmooth optimization techniques. 

3. Given an approximate or qptiiral rraximizer of (5.2), say ( u k+1 , , u N ) , 
let w 2 denote the optimal solution of the two-dimensional 
assignment problem (4.10) corresponding to this maximizer 

4. Formulate the recovery (N-k) -dimensional problem (4.19), 
modified as discussed at the end of subsection IV. 3. 3 for 
sparse problems . At this stage, z N , as defined in (4.21), 
contains the alignment of the indices {i lf ... , i^ +1 } . 

Formulate the final two-dimensional problem (4,22) and 
complete the final recovered solution z N as in (4.25). 

To explain the lower and upper bounds, let @ N be as 
defined in (4.8) with k=l f let V N {z N ) be the objective function 
value of the N- dimensional assignment problem in equation 
(3.5) corresponding to a feasible solution z N of the 
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constraints in (3.5), and let 2^~ k+1 be an optimal solution of 
(3.5). Then, 

0 N (u 3 ,...,u N ) < V N (z N ) z V N (2 N ) (5.3) 

is the desired inequality. 
COMMENTS : 

1. Step 2 is the computational part of the algorithm. 
Evaluating ^.^i and computing a subgradient used in the 
search procedure requires 99% of the computing time in 
the algorithm. This part uses two-dimensional assignment 
algorithms, a search over a large number of indices, and 
a nonsmooth optimization algorithm. It is the second 
part (the search) that consumes 99% of the computational 
time and this is almost entirely parallelizable . 

2. In track maintenance, the warm starts for the initial 
multipliers for the first step are available. For the 
relaxation procedure, initial multipliers are available 
in step 2 from the prior step in the algorithm. 

3. There are many variations on the above algorithm. For 
example, one can compute a good solution on the first 
pass {k-1) and not perform the optimization in (2) 
thereafter. This yields a great solution. Thus one can 



10 



CSUR.01USRI 

continue the optimization at the first pass, and 
immediately compute quality feasible solutions to the 
problem. 

VI. LAGRANGIAN RELAXATION ALGORITHM 
FOR THE THREE-DIMENSIONAL ALGORITHM 

In this section, we illustrate the relaxation algorithm 

for the three-dimensional assignment problem. Having 

discussed the general relaxation scheme, 

2J zJ C. . . Z. . . 

Vo vo vo W3 W3 
Subject To £ z i 1 ± 2 i 3 =1 f i 1 = l.»./^/ 

i 1 =0 i 3 =0 1 2 3 

. „ E Z W3 =1 ' i 3 = 1 — M 3/ 

i^O i 2 = 0 1 2 3 

Z W 3 ^ (0,1) ^or all 



(6.1) 



To ensure that a feasible solution of (6.1) always exists for 
a sparse problem, all variables with exactly one nonzero index 
(i.e., variables of the form z i 00 r z 0i0 r an d z QOi for i p =l,...,M p 
15 and p=l,2,3 are assumed free to be assigned with the 
corresponding cost coefficients being well-defined) . This 
assumption is valid in the tracking environment (A.B. Poore . 
Multidimensional assignments and multitarget tracking, 
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partitioning data sets, I.J. Cox, P. Hansen, and B. Julesz, 
editors, DIMACS Series in Discrete Mathematics and Theoretical 
Computer Science, American Mathematical Society, Providence, 
R.I., 19:169-198, 1995; A.B. Poore . Multidimensional 
5 assignment formulation of data association problems arising 
from multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications, 3:27-57, 1994). 



VI . 1 . The Lagrangian Relaxed Assignment Problem 
10 The three-dimensional assignment problem (6.1) has three 

sets of constraints and one can describe the relaxation by- 
relaxing any of the three sets of constraints, the description 
here is based on relaxing the last set of constraints. A 
(M 3 +l) -dimensional multiplier vector (i . e . , u 3 e ^ 3 +1 ) 
15 associated with the 3 rd constraint set will be denoted by 
u = lu 0 , u x , u M ) with u Q =0 being fixed for notational 
convenience. (The zero multiplier u 0 3 = 0 is used for notational 
convenience.) Now, the three-dimensional assignment problem 
(6.1) is relaxed to a two-dimensional assignment problem by 
2 0 incorporating last set of constraints into the objective 
function via the Lagrangian. Although any constraint set can 
be relaxed, we choose the last set of constraints for 
convenience. The relaxed problem is 



193 



CSUR.01USRI 



a, m, h, 



(6.2) 



<P 3 (u 3 ) = Minimize <p 3 (z 2 ; u 3 ) = £ £ £ c^^zf^ 

l t =0 i 2 = 0 i 3 =0 

+ E E E 4v 3 - 

i 2 =0 

Subject to 2l 2l z ± 1 i 2 i 3 = 1 f i 1 = l/-»/M L f 

i 2 =0 i 3 =0 



K. M, 

Minimize 22 z2 z2 

V° i 2 = 0 1 3 = 0 



1^3 f — % "3 

2 3 = 0 



i 1 =0 i 3 =0 



One of the major steps in the algorithm is the 
maximization of & 3 {u 2 ) with respect to the multiplier vector 
5 u 3 . It turns out that @ 3 is a concave, continuous, and 
piecewise affine function of the multiplier vector u 3 , so that 
the maximization of @ 3 is a problem of nonsmooth optimization. 
Since many of these algorithms require a function value and a 
subgradient of & 3 at any required multiplier value (u 3 ) , we 
10 address this problem in the next subsection. We note, 
however, that there are other ways to maximize <P 2 and the next 
subsection addresses but one such method. 

VI. 2. Properties Lagrangian Relaxed Assignment Problem 
15 For a function evaluation of <P 3/ we show that an optimal 

{or suboptimal) solution of this relaxed problem (6.2) can be 
constructed from that of a two-dimensional assignment problem. 
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Then, the nonsmooth characteristics of <P 3 are addressed, 
followed by a method for computing a subgradient. 

Evaluation of 0 3 . Define for each (i 1 ,i 2 ) an index 



(6.3; 



j 3 = j 3 (i lf i 2 ) and a new cost function c? la by: 

j 3 = j 3 (i lf i 2 ) =arg min{c* ii2i3 + ul | i 3 = 0, 1, ...,M 3 J, 
<i 2 -<i 2 j 3 (V 2 ) + ^3 for {i lf i 2 ) M0,0) ( 

c o 2 o = £ Minimum 0 , c 0 3 0i + u? . 

i 3 =0 1 3 3) 

Given an index pair (i 1 ,i 2 ) , j 3 need not be unique, 
resulting in the potential generation of several subgradients 
(6.11). (We further discuss this issue at the end of the 
Subsection. ) 
10 Now, define 

4(u 3 ) =Minimize <f 3 {z 2 ;u 3 ) ^£ £ <i 2 <i 2 ~£ 



Subject To $J z. 2 , =1, 1 =1,...,M,, 

Vo 12 1 1 

2 

^' ^2 1/ " " * '^2' 

z, 2 . e { 0 , 1 } for all L, i 0 . 

12 j- ^ 



(6.4) 



As an aside, two observations are in order. The first is that 
the search procedure needed for the computation of the relaxed 
cost coefficients in (6.3) is the most computationally 
intensive part of the entire relaxation algorithm. The second 
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is that a feasible solution z 3 of a sparse problem (6.1) yields 
a feasible solution z 2 of (6.4) via the construction 



2 

z i 1 i 2 



1/ if *l x i z ,i^ for some (i lf i 2 ,i 3 ), 
0, otherwise. 



5 Thus the two-dimensional assignment problem (6.4) has feasible 
solutions other than those with exactly one nonzero index if 
the original problem (6,1) does. 

The following Theorem 6.1 states that an optimal solution 
of (6.2) can be computed from that of (6.4) . Furthermore, if 
10 the solution of either of these two problems is e -optimal, 
then so is the other. The converse is contained in Theorem 
6.2. 

Theorem 6.1. Let w 2 be a feasible solution to the two- 
dimensional assignment problem (6.4) and define w 3 by 

15 3 2 (6.5) 

w V 2 i3 = \i 2 if i 3 = j 3 and i 2 ) * (0, 0) , 

w i x i 2 i 3 = 0 if and Ui^a* * (0 ' 0) ' 

*oo i 3 =0 if c 0 3 oi 3 + ^ 3 >0. 

Then w 3 is a feasible solution of the Lagrangian relaxed 
problem (6.2) and 0 ?{w 3 ; u 3 ) = <pAw 2 ; u 3 ) . If, in addition, w 2 is 
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optimal for the two-dimensional problem, then w 3 is an optimal 
solution of the relaxed problem and & 3 {u 3 ) =& 3 (u 3 ) . 

With the exception of one equality being converted to 
5 an inequality, the following theorem is a converse of 
Theorem 6.1. 

Theorem 6.2. Let w 3 be a feasible solution to problem (6.2) 
and define w 2 by 

3 * (6 - 6) 

S <i 2 = E <i 2 i 3 for (i lf i 2 ) * (0,0) , 

j| w£ 0 = 0 if (i lf i 2 ) = (0,0) and c^ i3 + u^>0 for all i 3 , 

* v Woo = 1 i f {i lf i 2 ) = (0,0) and c 0 3 oi3 + u^< 0 for some i 3 . 

W Then w 2 is a feasible solution of the problem (6.4) and 

0 3 (w 3 ;u 3 ) z <p 3 {w 2 ;u 3 ) . If, in addition, w 3 is optimal for (6. 2), 
then w 2 is an optimal solution of (6.4), <p 3 {w 3 ; u 3 ) = <f 3 {w 2 ; u 3 ) 
and 0 3 {u 3 ) =& 3 {u 3 ) . 

15 

An Algorithm for Computing @ 3 and a Subgradient. 

Given u 3 the problem of computing & 3 {u 3 ) and a subgradient 
of <P 3 {u 3 ) can be summarized as follows. 
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Starting with problem (6.1), form the relaxed problem 
(6.2) . 

To solve (6.2) optimally, define the two-dimensional 
assignment problem (6.4) via the transformation (6.3). 
Solve the two-dimensional problem (6.4) optimally. 
Reconstruct the optimal solution, say w 3 , of (6.2) via 
equation (6.5) as in Theorem 6.1. 
Then 



10 



0. 



(u 3 ) ^<p 3 (w 3 ;u 3 ) =}2 £ £ cl^iAi^ 
V° Vo V° 



i 1 =0 i 2 =0 



(6.7) 



A subgradient is given by substituting w into the 
objective function (6.2), erasing the minimum sign, and 
taking the partial with respect to u? . The result is 



9u 



d$ 3 (u 3 ) 

- L -i 



i,=0 i,=0 1 2 3 



for 1=1, 



(6.8) 



VI . 3 . The Recovery Procedure 
15 The next objective is to explain a recovery procedure, 

i.e., given a feasible (optimal or suboptimal) solution w 2 of 
(6.4) (or w 3 of (6.2) constructed in Theorem 6.1), 

generate a feasible solution z 3 of equation (6.1) which is 
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close to w 2 in a sense to be specified. We first assume 
that no variables in (6. 1) are preassigned to zero; this 
assumption will be removed shortly. The difficulty with the 
solution w 3 in Theorem 6.1 is that it need not satisfy the 
5 third set of constraints in (6.1). (Note, however, that if w 2 
is an optimal solution for (6.4) and w 3 , as constructed in 
Theorem 6.1, satisfies the relaxed constraints, then w 3 is 
optimal for (6.1).) The recovery procedure described here is 
designed to preserve the 0-1 character of the solution w 2 of 

2 

10 (6.4) as far as possible: If ^1=1 and i^O or i 2 *0, the 
corresponding feasible solution z 3 of (6.1) is constructed so 
that z ii± =l for some (i 1 ,i 2 ,i 3 ). By this reasoning, variables 
of the form z$ 0i ^ can be assigned a value of one in the 
recovery problem only if w£ 0 = l . However, variables Zq 01 ^ will 

15 be treated differently in the recovery procedure in that 
they can be assigned either zero or one independent of the 

2 

value w 00 . This increases the feasible set of the recovery 
problem, leading to a potentially better solution. 

Let |(i a (1 2 ) 7 i 2 (I 2 )|^ ^ be an enumeration of indices {i 1 ,i 2 } 
2 0 of w 2 (or the first 2 indices of w 3 constructed in equation 
(6.5) ) such that w^^^^ =1 for (i x (1 2 ) , i 2 (1 2 ))* (0, 0) and 
(i 1 (1 2 ) , i 2 (1 2 ) )= (0, 0) for I 2 =0 regardless of whether w 0 2 0 = l or 
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not. (The latter can only improve the quality of the feasible 
solution. ) 

Next to define the two-dimensional assignment problem 

that restores feasibility, let 

2 3 (6.9) 
Ci 2 i 3 = c ii(l2)i2(l2)i3 for l^O,...,^. 

Then the two-dimensional assignment problem that restores 
feasibility is 



v 2 ^ 2 2 



(6.10) 



Minimize 

Vo "o ^ 3 ^ 3 



2 

Subject to 2^ z 2 i =1, I 2 = 1,...,L 2 , 



2-3 - 3 " 3 

z, 2 , d0,l} Tor all I 9 ,i 



2' ^3' 



10 The next objective is to remove the assumption that all 

cost coefficients are defined and all zero-one variables are 
free to be assigned. We first note that the above 
construction of a reduced two-dimensional assignment problem 
(6.11) is valid as long as all cost coefficients c 2 are defined 

15 and all zero-one variables in z 2 are free to be assigned. 
Modifications are necessary for sparse problems. If the cost 
coefficient c^a^ ± 2 (i 2 ) ± 3 is well-defined and the zero-one 
variable Zi^i^i^i^^ is free to be assigned to zero or one, then 
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define c\^= c?^,^,^ as in (6.9) with zj^,^,^ being free 
to be assigned. Otherwise, z 1± is preassigned to zero or the 
corresponding arc is not allowed in the feasible set of arcs. 
To ensure that a feasible solution exists, we now only need 
ensure that the variables z^ Q are free to be assigned for 
12 = 0,1,...,!^ with the corresponding cost coefficients being 

3 3 

well-defined. Recall that all variables of the form z^ 00 , z^ t ^ 



and z O0i3 are free (to be assigned) and all coefficients of the 



3 3 3 

corresponding form c i oo , c oi ^ 0 and c 00i ^ need to be defined. This 



10 is accomplished as follows: if the cost coefficient c i (I }i (i )0 



is well-defined and z i (i )i {i )o ^ s free, then define 

2 3 2 

°i 2 o = c i 2 a 2 )i 2 (i 2 )o with z 1Q being free. Otherwise, since all 
variables of the form z? 0 and Zq 10 are free to be assigned 
with the corresponding costs being well-defined, set 

2 3 3 

15 c i 2 o = c i 1 (i 2 )oo + c oi 2 (i 2 )o / where the first term is omitted if i 1 (l 2 )=0 
and the second, if i 2 (2 2 )=0. For i 1 (l 2 )=0 and i 2 (2 2 )=0 , define 

2 _ 3 

c oo ~ c ooo • 

VI. 4. The final recovery problem 
The recovery problem for the case N = 3 is 
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(6.11) 



Subject to 2J z ii =1/ ^ 2 = -^2' 
i 3 =o 2 3 

1 2 =0 2 3 

z i 2 i 3 e {0,1} for all I 2 ,i 3 . 

Let Y be an optimal or feasible solution to this two- 
dimensional assignment problem. The recovered feasible 
5 solution z 2 is defined by 

^3 fl, if WI^Wl,) and 3^=1, (6 ' 12) 
2li2i3 lo, otherwise. 

Said in another way, let {(1 2 (1 3 ) , i 3 (2 3 )}^ 3 =0 be an enumeration 
of indices {l 2 ,i 3 } of Y such that if . =1 for 
10 (1 2 (1 3 ) ,i 3 (l 3 ))*(0,0) and (I 2 (1 3 ) , i 3 (1 3 ) ) = (0 , 0) for I 3 =0 
regardless of whether Yoo =1 or not. Then the recovered 
feasible solution can be written as 



^3 fl/ if i 1 = i 1 (I 12 ) , i 2 = i 2 (i 12 ) ,i 3 = l 3 (I 3 ) , 
i i 1 2 i 3 [o, otherwise. 

15 

The upper bound on the feasible solution is given by 

M, M~ M, 



^ 3 ) = £ £ £ c wA 2 i 3 (6 - 14) 

i^O i 2 = 0 i 3 = 0 123 123 

20 and the lower by €> 3 {u 3 ) for any multiplier value (u 3 ) . 
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In particular, we have 

0 3 (u 3 ) <$ 3 (u 3 ) <; V 3 (J 3 ) tV 3 (z 3 ) , 
where u 3 is any multiplier value, if is a maximizer of @ 3 {u 3 ) , 
—3 is an optimal solution of the multidimensional assignment 
5 problem and z 3 is the recovered feasible solution defined by 
(6.13) . 

VII. OTHER RELAXATIONS FOR THE MULTIDIMENSIONAL 
ASSIGNMENT PROBLEM 

In this section, we briefly describe other possible 

10 relaxations and their implications. These include linear 

programming relaxations and the corresponding duals. 

Recall from Section II that one starts with the 

definition of the zero-one variable z. . =1 and cost 

1 n" :l n 

coefficients to form the problem 

15 



Minimize c. . z. . 

- L i'*- L w 

Subject To X) z i...i =1 (^ = 1,-,^), 

7 7 -7' IN 

(i 2 = l,~,M 2 ), (7.1) 

1* 7 ... T IN 



. .2, . Vv 1 

2 0 ifViVr 1 * 



(i p = l,-,M p and p=2,~,W-l), 



v/f 0 ' 1 ! for a11 
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where c 0 0 is arbitrarily defined to be zero. Here / each group 
of sums in the constraints represents the fact that each non- 
dummy report occurs exactly once in a "track of data" . One 
can modify this formulation to include multiassignments of 
5 one, some, or all the actual reports. The assignment problem 
(3.5) is changed accordingly. For example, if z* is to be 
assigned no more than, exactly, or no less than n± times, then 

k 

the w = l" in the constraint (3.5) is changed to "£,=,^12*," 
respectively. Modifications for group tracking and 

10 multiresolution features of the surveillance region will be 
addressed in future work. In making these changes, one must 
pay careful attention to the independence assumptions that 
need not be valid in many applications. 

Again, the recent work of Poore and Drummond (A.B. Poore 

15 and O.E. Drummond. Track initiation and maintenance using 
multidimentional assignment problems. In D.W. Hearn, W.W. 
Hager, and P.M. Pardalos, editors, Network Optimization, 
volume 450, pages 407-422, Boston, 1996. Kluwer Academic 
Publishers B.V.) significantly extends the formulation of the 

20 track maintenance and initiation to new approaches. The 
discussions of this section apply equally to those formu- 
lations . 

VI 1 . 1 . The Linear Programming Relaxation 
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In the linear programming relaxation, one replaces the 
zero-one variables constraints 

z. 6{0,1} for all i lf ...,i N (7.2) 

with the constraint 
5 0<z. . <1 for all i 1f .„i„. (7.3) 

Then, the problem (7.1) can be formulated as a linear 
programming problem with the constraint (7.2) in (7.1) 
replaced by (7.3) with a special block structure to which the 
Dant zing-Wolf e decomposition is applicable. Of course, after 
10 solving this problem, one must now recover the zero-one 
character of the variables in (7.1) and there are many ways to 
do this, such as using the two-dimensional assignment 
problems. Commercial software is also available. 

15 VII. 2. The Linear Programming Dual and 

Partial Lagrangian Relaxations 

Given the linear programming relaxation, one can formulate 

the dual problem or the partial Lagrangian relaxation duals 

with respect to any number of constraints. In particular, 

20 this is precisely what is done in Section 3 on the Lagrangian 

relaxation algorithm presented. The much broader class of 

algorithms provided in the US Patent 5537119 (Aubrey B. Poore, 

Jr., Method and System for Tracking Multiple Regional Objects 

by Multi-Dimensional Relaxation, US Patent Number 5537119, 
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filed 14 March 1995, issue date of 16 July 1996) can be 
modified to remove the zero-one character when one relaxes M 
sets of constraints to an (N-M) -dimensional problem and 
recovers with an (N-M +1) -dimensional problem. This avoids 
5 the problems associated with the NP-hard (N-M) - and (N-M +1) - 
dimensional problems. However, to restore the zero-one 
character, one can do it sequentially with an assignment 
problem or with one of the many zero-one rounding techniques. 
These formulations are easy to work out and thus the details 

10 are omitted. 

The complete dual problem is another way of solving the 
LP problem and may indeed be more efficient in certain 
applications. In addition, the solution of this dual problem 
may provide excellent initial approximations to the 

15 multipliers for Lagrangian relaxations. 



VIII. HOT STARTS FOR TRACK MAINTENANCE 

AND INITIATION: BUNDLE MODIFICATIONS 

Thus reuse of multipliers and the first proof that this 

2 0 reuse is actually valid was presented in this section. The 

reuse in track maintenance is demonstrated and discussed in 

the work of Poore and Drummond (A.B. Poore and O.E. Drummond. 

Track initiation and maintenance using multidimensional 

assignment problems. In D.W. Hearn, W.W. Hager, and P.M. 
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Pardalos, editors, Network Optimization, volume 450, pages 
407-422, Boston, 1996. Kluwer Academic Publishers B.V.). The 
only item left is the modification of the bundle of 
subgradients for the use with the multiplier values as one 
5 goes through the recovery problem as in Section IV. 3. 6 or as 
one moves the window in track maintenance. It is this aspect 
of the nonsmooth optimization that adds an order of magnitude 
to the speed of the algorithms. This work is based on both a 
mathematical proof as in Section IV. 3. 6 as well as 

10 computational experience and heuristics. 

IX. ADAPTIVE AND VARIALBE WINDOW SIZE SELECTION 
These data association algorithms, based on the 
multidimensional assignment problem, range from sequential 
processing to many frames of data processed all at once! The 

15 data association problem for multisensor-multitarget tracking 
is formulated using a window sliding over the frames of data 
from the sensors. This discussion focuses on the work of 
Poore and Drummond and not on the formulations in Section III. 
Firm data association decisions for existing track are made 

2 0 at, say frame k, with the most recent frame being k+M. 
Decisions after frame k are soft decisions. Reports not in 
the confirmed tracks are used to initiate tracks over frames 
numbered k-I to k+M. 
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The length of these windows varies from sequential 
processing, which corresponds precisely to J = 0 and M = 1, to 
user defined large values of I and M. There is a marked change 
in performance over this range. As the two windows become 
5 exceedingly long, there is little statistically significant 
information gained and indeed performance can degenerate 
slightly. Thus, one must manually find the correct window 
lengths for performance in a given scenario and the algorithms 
do not change for a changing environment . Thus we propose an' 
10 adaptive method for adjusting the window lengths. (The method 
has been highly successful for selecting the correct window 
length for multistep methods in ordinary differential 
equations . ) 

1. Compute the solution for different window lengths that 
15 overlap and differ by one or two frames. 

2. Compare the solution quality, i.e., the value of the 
objective function, for two adjacent window lengths. 

20 3. If there is a predefined improvement in either direction, 
we then, for stability, repeat the exercise for a shorter 
or longer on the first try. If there is consistent 
information, we adjust the window size in the indicated 
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direction. This can be given a well defined mathematical 
formula in terms of the assignment problems of different 
dimensions . 



X. SENSITIVITY ANALYSIS 
For sequential processing in which the two-dimensional 
assignment algorithm is used, one can use the LP sensitivity 
analysis theory and easily obtain the corresponding answers. 
For the multiframe processing, the optimal solution is not 
available; however, there are still two approaches one can 
use. First, the basic algorithm can perform the same 
sensitivity analysis at each stage (loop) of the algorithm as 
is done in the two-dimensional assignment problem since the 
evaluation of function $ N _ k+1 is equivalent to a two-dimensional 
assignment problem. Alternately, one could use an LP 
relaxation of the assignment problem and base the sensitivity 
on the resulting LP problem. We currently see this as an 
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important step in finding even better solutions to the 
assignment problem if so desired. 



XI . NEW AUCTION ALGORITHMS . 
In this section we present a new auction algorithm for the 
two-dimensional assignment problem. 

An important step in solving N-dimensional assignment 
problems for N > 3 is finding the optimal solution of the two- 
dimensional assignment problem. In particular, we wish to 
solve the two-dimensional problem which includes the zero- 
index. Typically this problem can be thought of as trying to 
find either a minimum cost or maximum utility in assigning 
person to objects, tenants to apartments, or teachers to 
classrooms. We will follow the work of Bertsekas and call the 
first index set persons and the second index set objects. The 
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statement of this problem is given below (11.1) when the 
number of persons is m and the number of objects is n. 



Minimize 7,7, ex.. 

*7-i *7—t 11 11 

1=0 j=0 J J 



Subject to £ x„ = l for all i = l,...,m, 

J=0 



(11.1) 



S x n = 1 for all j = l, . . . ,n. 

i=0 

There are a couple of assumptions which we will make 
about (11.1). First of all, we assume that c i0 and c oj are 
well-defined and the corresponding variables x i0 and x oj are 
free to be assigned. Second if a cost c ±j happens to be 
undefined, then the corresponding variable x ±j will be set to 
0. In effect if c ±j is undefined, we simply remove this cost 
and variable from inclusion in the problem. 

Notice that there are no constraints on the number of 
times person 0 and object 0 can be assigned. But notice that 
the first constraint requires that each non-zero person i must 
be assigned exactly once. Similarly, the second constraint 
forces each non-zero object j to be assigned exactly one time. 
Finally, the last constraint gives an integer solution, 
although we will see shortly that this constraint can be 
relaxed to admit any solution x ±j *0 . One reason for requiring 
that all of the costs of the form c i0 and c oj be defined is so 
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that we are guaranteed a feasible solution exists for the 
given problem. 

XI. 1. Relaxation of One Constraint 
Relaxing the last constraint, we obtain the following 
problem: 



Min 



inize £ £ c x + £ u] 

i = 0 j = 0 J J j=0 



2=0 j=0 

i=0 j=0 



C H* U 3 



2 = 0 ^ 



j = 0 



Subject to £ x^.-l for all 2 = 1,..., 



j=o 

X^6{0,1}. 

This lower bound is achieved and the second constraint in 
the original problem is satisfied if the following conditions 
are met : 



For i=0, 



For i*0, 



x oj = 



1, if c 0j + u]<0\ 

0, otherwise 

1, if j = argmin{c ik +u£\k=0,l,...,n }\ 
0 otherwise 



All of j f s are assigned only one time. 

The relaxation of the first constraint is analogous and would 
lead to similar results with i and j exchanged and the 
introduction of the multiplier u 1 . 

XI. 2. Relaxation of Both Constraints 
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This time we will relax both constraints and use duality 
theory to obtain a relationship between the multipliers u 1 and 
u 2 . 

Definition: An assignment S and a multiplier set (u 1 ,u 2 ) are 
said to satisfy 6 -complementary slackness (E - CS) if 

c ± . + ul + Uj=0 for all (i,j)eS, 
c ij + u i +u j 2 ^" e for a11 (i,j)€A. 

Forward Auction Algorithm 

(1) Select any unassigned person 1 

(2) Determine the following quantities: 

j ± = arg min j c lk +u£\k eA(i)}, 

w.= min | c. k +ul\k EA[i) , k*j^. 

In the selection of j ± above, if a tie occurs between 0 
and any non-zero index k, then select j x as k. Otherwise, 
if there is a tie between two or more non-zero indices, 
the choice of j 1 is arbitrary. 

Also if A(i) consists of only one element, then set w ± =°° . 

(3) Update the multipliers and the assignment: 
If Ji = 0, then 
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(a) Add (1,0) to the assignment. 

(b) Update u a x : = -c icr 
If j x # 0, then 

(a) Add (i/Ji) to the assignment. 

(b) Remove from the assignment if j t was 
previously assigned. 



(c) Update u 7 2 := + (w.-v.) +e = w.-c. +e. 

Ji Ji 2 2' 2 2J i 

(d) Update u 2 2 := -^o ± . , + = -w\-e. 



Reverse Auction Algorithm. 

(1) Select any unassigned object j, such that c Qj + Uj >0, 

(2) Determine the following quantities: 



i^arg min | c JJt + Ujf | Jr eB(j)}, 



In the selection of above, if a tie occurs between 0 
and any non-zero index k, then select j ± as k. Otherwise, 
if there is a tie between two or more non-zero indices, 
the choice of j i is arbitrary. 
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Also if B{j) consists of only one element, then set 
Yj = 00 - 

(3) Update the multipliers and the assignment: 
If ij = 0, then 

(a) Add (0,j) to the assignment. 

(b) Update Uji = -c Q .. 
If ij * 0, then 



(a) Add (ij,j) to the assignment. 

(b) Remove (i-^jO from the assignment if (i jf j) was 
previously assigned. 

(c) Update 12^1 = 11^ + (Yj-Pj) + e = Y-, - ^ j + e. 

(d) Update u j 2 : = -( Cl + u^)= -y^-e . 
Combined Forward/Reverse Auction Algorithm 

1. Assume that u 2 is given as an arbitrary multiplier. 

2. Adjust the value of u 2 for each object j as follows: 

If c Qj + u] <0, then set u* : = -c Q _. . 

3 . Run iterations of the Forward Auction Algorithm 
until all persons become assigned. 

4. Run iterations of the Reverse Auction Algorithm 
until all of the objects become assigned. 

Note at the completion of the Forward auction step we 
have the following conditions satisfied: 

215 



CSUR.01USRI 

• c 0J + u? £ 0 for all objects j. 

• c^. + u^ + u? =0 for all (i,j)€ S. 

• c ij + u j 2 * m i n { c ijt +u Jt 2 } + e for all ( i , j ) £ S. 
Thus we can prove the following proposition. 

2 

5 Proposition : If we assume that c Q _. + ^0 at the start of the 
Forward Auction Algorithm and all of the persons are assigned 
via a forward step, then we have: 

c^ + ul + Uj > -E for all A. 

c.. + ul + uj =0 for all S. 

10 c ±j + u j * min | c ijc +u M +e for all S. 



Optimality of the Algorithm 

Theorem: 6-CS preserved during every forward and reverse 
iteration. 

15 Theorem: If a feasible solution exists, then the resulting 
solution is with m€ of being optimal for the 
Combined Forward Reverse Algorithm. 



XII. SOME CONCLUDING COMMENTS 
2 0 Although the algorithm appears to be serial in nature, 

its primary computational requirements are almost entirely 
parallelizable . Thus parallelization is planned. 
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Step 2 is the computational part of the algorithm. 
Evaluating ® N - k+1 and computing a subgradient use in the search 
procedure requires 99% of the computing time in the algorithm. 
This part uses two-dimensional assignment algorithms, a search 
5 over a large number of indices, and a nonsmooth optimization 
algorithm. It is the second part (the search) that consumes 
99% of the computational time and this is almost entirely 
parallelizable . Indeed, there are two-dimensional assignment 
solvers that are highly parallelizable. Thus, we need but 

10 parallelize the nonsmooth optimization solver to have a 
reasonably complete parallelization. 

If a sensitivity analysis is desired or if one is 
interested in computing several near-optimal solutions, a 
parallel processor with a few powerful processors and good 

15 communication such as on the Intel Paragon would be most 
beneficial . 

The foregoing discussion of the invention has been 
presented for purposes of illustration and description. 
Further, the description is not intended to limit the 
2 0 invention to the form disclosed herein. Consequently, 
variation and modification commensurate with the above 
teachings, within the skill and knowledge of the relevant art, 
are within the scope of the present invention. The embodiment 
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described hereinabove is further intended to explain the best 
mode presently known of practicing the invention and to enable 
others skilled in the art to utilize the invention as such, or 
in other embodiments, and with the various modifications 
required by their particular application or use of the 
invention. It is intended that the appended claims be 
construed to include alternative embodiments of the invention 
to the extent permitted by the prior art . 
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